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Mastering Differential Equations: The Visual Method 


Scope: 

T he field of differential equations goes back to the time of Newton 
and Leibniz, who invented calculus because they realized that 
many of the laws of nature are governed by what we now call 
differential equations. 

A differential equation is an equation involving velocities or rates of change. 
More precisely, it’s an equation for a missing mathematical expression or 
expressions in terms of the derivatives (i.e., the rates of change) of these 
expressions. Differential equations arise in all areas of science, engineering, 
and even the social sciences. The motion of the planets in astronomy, the 
growth and decline of various populations in biology, the prediction of 
weather in meteorology, the back-and-forth swings of a pendulum from 
physics, and the evolution of chemical reactions are all examples of 
processes governed by differential equations. 

In the old days, when we were confronted with a differential equation, the 
only technique available to us to solve the equation was to find a specific 
formula for the solution of the equation. Unfortunately, that can rarely be 
done—most differential equations have no solutions that can be explicitly 
written down. So in the past, scientists would often come up with a simpler 
model for the process they were trying to study—a model that they could 
then solve explicitly. Of course, this model would not be a completely 
accurate description of the actual physical process, so the solution would 
only be valid in a limited setting. 

Nowadays, with computers (and, more importantly, computer graphics) 
readily available, everything has changed. While computers still cannot 
explicitly solve most differential equations, they can often produce 
excellent approximations to the exact solution. More importantly, computers 
can display these solutions in a variety of different ways that allow 
scientists or engineers to get a good handle on what is happening in the 
corresponding system. 
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Scope 


In this course, we take this more modem approach to understanding 
differential equations. Yes, we grind out the explicit solution when possible, 
but we also look at these solutions geometrically, plotting all sorts of 
different graphs that explain the ultimate behavior of the solutions. 

In each of the 4 sections of this course, we begin with a simple model and 
then use that model to introduce various associated topics and applications. 
No special knowledge of the field from which the model arises is necessary 
to understand the material. Two broad subthemes appear over and over 
again in the course: One is the concept of bifurcation—how the solutions 
of the equations sometimes change dramatically when certain parameters 
are tweaked just a little bit. The second subtheme is chaos—we will see 
how extremely unpredictable behavior of solutions occasionally arises 
in the setting of what would seem to be a relatively simple system of 
differential equations. 

In the first part of the course, we discuss the simplest types of differential 
equations—first-order differential equations—using several population 
models from biology as our examples. The simplest model is the unlimited 
population growth differential equation. Later examples, such as the logistic 
population growth model (a limited growth model), take into account the 
possibility of overcrowding and harvesting. Using elementary techniques 
from calculus, we solve these equations analytically when possible and plot 
the corresponding qualitative pictures of solutions, such as the slope field, the 
phase line, and the bifurcation diagram. We also describe a simple algorithm 
that the computer uses to generate these solutions and show how this 
technique can sometimes fail (often because of chaos). Finally, we introduce 
a parameter into our model and see how very interesting bifurcations arise. 

In the second part of the course, we turn our attention to second-order 
differential equations. Our model here is the mass-spring system from 
physics. We see relatively straightforward behavior of solutions as long as 
the mass-spring system is not forced, but when we introduce periodic forcing 
into the picture, the much more complicated (and sometimes disastrous) 
behaviors known as beating modes and resonance occur. 
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The next part of the course deals with systems of differential equations. Our 
model here is the predator-prey model from ecology. At first we concentrate 
on linear systems. 

We see how various ideas from linear algebra allow us to solve and analyze 
these types of differential equations. Bifurcations recur as we investigate 
the trace determinant plane. We then move on to a topic of considerable 
interest nowadays, nonlinear systems. Almost all models that arise in nature 
are nonlinear, but these are the differential equations that can rarely be 
solved explicitly. We investigate several different models here, including 
competitive systems from biology and oscillating chemical reactions. 
Another model will be the Lorenz system from meteorology; back in the 
1960s, this was the first example of a system of differential equations that 
was shown to exhibit chaotic behavior. 

In the final part of the course, we turn to the concept of iterated functions 
(also called difference equations) to investigate the chaos we observed in 
the Lorenz system. Our model here is an iterated function for the logistic 
population model from biology, a very different kind of model than our 
earlier logistic differential equation. In this case, we see that lots of chaos 
emerges, even in the simple iterated function, and we begin to understand 
how we can analyze and comprehend this chaotic behavior. 

While calculus is a central notion in differential equations, we will not delve 
into many of the specialized techniques from calculus that can be used to 
solve certain differential equations. The only topics from calculus that we 
presume familiarity with are the derivative, the integral, and the notion of a 
vector field. Any other relevant concepts from calculus or linear algebra are 
introduced before being used in the course. 

By the end of this course, you will come to see how all the concepts from 
algebra, trigonometry, and calculus come together to provide a beautiful and 
comprehensive tool for investigating systems that arise in all areas of science 
and engineering. You will also see how the field of differential equations is 
an area of mathematics that, unlike algebra, trigonometry, and calculus, is 
still developing, with many new and exciting ideas sparking interest across 
all disciplines. 
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Lecture 1: What Is a Differential Equation? 


What Is a Differential Equation? 

Lecture 1 


A n ordinary differential equation (ODE) is an equation for a 
missing function (or functions) in terms of the derivatives of those 
functions. Recall that the derivative of a function y(t) is denoted 
either by y'(t) or by dyldt. Suppose y(t 0 ) =yo . Then the derivative // 0 ) gives 
the slope of the tangent line to the graph of the function y(t) at the 
point (t 0 ,y 0 ). 

In the old days, the only tools we had to solve ODEs were analytical 
methods—a variety of different tricks from calculus that sometimes enabled 
us to write down an explicit formula for the solution of the differential 
equation. Unfortunately, most ODEs cannot be solved in this fashion. But 
times have changed: Now we can use the computer to approximate 
solutions of ODEs. And we can use computer graphics and other geometric 
methods to display solutions graphically. So this gives us 2 new ways to 
solve differential equations, and these are the methods that we will 
emphasize in this course. 

Probably the best-known differential equation (and essentially the first 
example of a differential equation) is Newton’s second law of motion. Drop 
an object from your rooftop. Ify measures the position of the center of mass 
of the object, then we would like to know its position at time t, that is, y(t). 
Newton’s law tells us that mass times acceleration is equal to the force on 
the object. So, if m is the mass, then we have my” = F(y ), where F is the 
force acting on the object when it is in position y(t). So we have a 
differential equation fory(f). 

Another example of a differential equation is the mass-spring system (or 
harmonic oscillator). Here y(t) is the missing function that measures the position 
of a mass attached to a spring attached to a ceiling. When we pull the mass 
down and let it go, y(t) gives us the resulting motion of the mass over time. The 
function y(t) is determined by the differential equation/' + by'+ ky= G(t ). 
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In this case, we have parameters like b and k as well as a forcing term G(t ) 
that depends only on t. We will see very different behaviors for the mass¬ 
spring system depending on these parameters and the forcing term. 
Sometimes a small change in one of these parameters creates a major 
change in the behavior of solutions. Such phenomena are called 
bifurcations, one of the subthemes we will encounter often in this course. 

Most of our course will be spent considering systems of differential 
equations. Systems of ODEs involve more than one missing function in the 
equation. We will consider numerous such examples, but the most famous 
is undoubtedly the Lorenz system from meteorology, below. 

x' = -10x + 10y 
y ' = -xz + Rx - v 
z' = xy-%z 

The Lorenz system of equations was the first example of a system of ODEs 
that exhibits chaotic behavior, another subtheme that we will encounter 
throughout the course. 

Let’s look at our first differential equation: the unlimited population growth 
model from biology, which is perhaps the simplest nontrivial differential 
equation. Suppose we have a species living in isolation (with no predators, 
no overcrowding, and no emigration) and we want to predict its population 
as a function of time. Call the population y(t). 

Our assumption is that the rate of growth of the population is directly 
proportional to the current population. This translates to the ODE y' = ky. 
Here k is a constant (a parameter) that depends on which species we are 
considering. This is an example of a first-order ODE, since the equation 
depends on at most the first derivative ofy(^). Usually we wish to find the 
solution of an initial value problem, that is, a specific solution of the ODE 
that satisfies y(0) = y 0 where y 0 is the given initial population. 

For simplicity, suppose k = 1, so our differential equation is y f = y. One 
solution is the exponential function y(t) = e\ since the derivative of the 
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Lecture 1: What Is a Differential Equation? 


function e is e l (so that y' does indeed equal y). Another solution is 
y(t) = Ce\ where C is some constant. Note that when C = 0, the solution is 
the constant function y(t) identically equal to zero. This is an equilibrium 
solution, one of the most important types of solutions of differential 
equations. The general solution of the equation is y(t) = Ce\ since we can 
solve any initial value problem using this formula. That is, given y( 0) = y 0 , 
the solution satisfying this initial condition is given by setting C = y 0 . 

More generally, we can consider the equation y f = ky, where k is some constant, 
say k = 2. As before, other solutions are of the form Ce 2t , where C is an arbitrary 
constant. Again we see that the solution with C = 0 is an equilibrium solution. 
We can solve any initial value problem by choosing C as our initial value y(0), 
so Ce 2t is the general solution of this differential equation. 

Here are 2 of the simplest methods for visualizing solutions of 
differential equations. 

1. The right-hand side of the differential equation tells us the slope of the 
solution of the ODE at any time t and population y. So in the t-y plane, 
we draw a tiny straight line with slope equal to the value on the right- 
hand side. A collection of such slopes is the slope field below. 


Figure 1.1 
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2. Then a solution must be everywhere tangent to the slope field, so we 
can sketch in the graphs of our solutions. 


Figure 1.2 



Note the constant solution, y(t) = 0. This is our equilibrium solution. 


Important Terms 


bifurcation: A major change in the behavior of a solution of a differential 
equation caused by a small change in the equation itself or in the parameters 
that control the equation. Just tweaking the system a little bit causes a major 
change in what occurs. 

equilibrium solution: A constant solution of a differential equation. 

general solution: A collection of solutions of a differential equation from 
which one can then generate a solution to any given initial condition. 

initial value problem: A differential equation with a collection of special 
values for the missing function such as its initial position or initial velocity. 
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Lecture 1: What Is a Differential Equation? 


ordinary differential equation (ODE): A differential equation that 
depends on the derivatives of the missing functions. If the equation depends 
on the partial derivatives of the missing functions, then that is a partial 
differential equation (PDE). 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap 1.1. 
Guckenheimer and Holmes, Nonlinear Oscillations. 

Hirsch, Smale, and Devaney, Differential Equations , chap 1.1. 
Roberts, Ordinary Differential Equations, chap 1.1-1.3. 
Strogatz, Nonlinear Dynamics and Chaos , chap 1.1. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , HPG Solver. 


Problems 


1. Let’s review some ideas from calculus that we used in this lecture. 

a. Compute the derivative of y(t) = t 3 + e\ 

b. Compute the second derivative of y(t) = t 3 + e*. 

c. Find a function^) whose derivative is t 3 + e t . 









d. Sketch the graph of the exponential function y(t) = e. 

e. Find the solution of the equation e (2?) = 1. 


2 . 


Sketch the slope field and solution graphs for the differential equation 
/=!• 


3. What are some solutions of the differential equation in problem 2? 


4. Repeat problems 2 and 3 for y' = t. 


5. In the unlimited population growth model y' = ky, what happens to 
solutions if k is negative? 


6. Find the general solution of the differential equation;;' = t. 


7. Find all of the equilibrium solutions of the differential equation 
y' = y 2 - 1. Also plot the slope field. What do you think will happen to 
solutions that are not equilibria? 


8. What would be the general solution of the simple differential equation 
y' = 0? What is the behavior of all of these solutions? 


9 


Lecture 1: What Is a Differential Equation? 


9. Consider the differential equation given by Newton’s second law, 
my" = g, where we assume that m and g are constants. Can you find 
some solutions of this equation? 


Exploration 


Think about the fact that back in the 1600s, Isaac Newton was able to 
come up with not only the second law of motion in physics but also the 
basics of calculus and differential equations (plus everything else he 
discovered in the sciences). He was quite an intellectual! There is plenty to 
read about his life and discoveries on the web. For an introduction, go to 
http://www.newton.ac.uk/newtlife.html. 
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A Limited-Growth Population Model 

Lecture 2 


W e begin this lecture by investigating a more complicated (but 
more realistic) population model, the limited-growth population 
model (also known as the logistic population growth model). The 
corresponding slope field gives us an idea of how solutions will behave, at 
least from a qualitative point of view. 

This does not always happen, however. Look at the slope field for the 
differential equation y' = y 2 - t, and put a little target at (2, 2). Can you see 
where a solution in the lower left quadrant should begin so that it hits the 
given target? 


Figure 2.1 
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Lecture 2: A Limited-Growth Population Model 


Well, maybe not. Nearby solutions tend to veer off from one another as y 
and t increase, as can be seen from the slope field, so the slope field does 
not tell us everything in this case. 


Figure 2.2 


y 



For the limited population growth model, we will assume that overcrowding 
may occur, which will hinder population growth or lead to population 
decline (if the population is too large). We make 2 assumptions about our 
population, y(t)\ 


• If y(t) is small, the rate of population growth is proportional to the 
current population (as in our unlimited growth model). 
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• There is a carrying capacity N (an “ideal” population size) such that 


y(t) > N means thaty(^) decreases ( y\t ) < 0), and 

y(t) < N means that y(t) increases ( y'(t ) > 0). 

The first assumption tells us that y'(t) = ky if y(t) is small. So we let our 
differential equation assume the form 




where the expression (?) 

1. is approximately equal to 1 if y is small; 

2. is negative if y > N\ 

3. is positive if y<N; and 

4. is 0 if y = N. 

Then one possibility (the simplest) for the term (?) is (1 —yIN). 

This yields the limited population growth ODE 

~ =ky( 1 -y/N). 
at 

For simplicity, let’s consider the case where k = N= 1: 

This doesn’t mean our carrying capacity population is just one individual 
(that would not be a very interesting environment). Rather, think of y(t) as 
measuring the percentage of the carrying capacity. We will solve this 
differential equation analytically in a later lecture, but for now we will use 
qualitative methods to view the solutions. 
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Lecture 2: A Limited-Growth Population Model 


Here is the slope field for the limited population growth model equation. 
Note the horizontal slopes when y = 0 or y « 1; these are our 
equilibrium points. 


Figure 2.3 



Here are some solution graphs. Note that they do exactly as we expected; 
they all tend to the equilibrium y = 1, the carrying capacity (assuming, of 
course, that the population is nonzero to start). 
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Figure 2.4 


P 



With an eye toward what comes later—when we deal with differential 
equations whose solutions live in higher dimensional spaces—we will plot 
the phase line for this ODE. This is a picture of the motion of a particle 
along a straight line where the position of the particle at time t is the value 
ofv(0- 
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Lecture 2: A Limited-Growth Population Model 


Figure 2.5 


W 

1 o 


/V 
o i» 


This course deals primarily with autonomous differential equations. These 
are ODEs of the form / = F(y) (i.e., the right-hand side does not depend on 
/). We simply find all the equilibrium points by solving F(y) = 0. Between 2 
such equilibria, the slopes are either always positive or always negative 
(assuming F(y) is a continuous function), so solutions either always increase 
or always decrease between the equilibrium solutions. For example, for the 
ODE / = y 3 - y, we have 3 equilibria at y = 0, 1, and -1. Above y = 1, we 
have / > 0. Between y = 0 and y = 1, we have / < 0. Between y = — 1 and 
y = 0, we have again / > 0. And below y = — 1, we have / < 0. So we know 
what the solutions will do, at least qualitatively. The phase line and some 
representative solution graphs are shown below. 


Figure 2.6 
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We can use the graph of the right-hand side of the differential equation to 
read this behavior off. Consider / = y 2 - 1. Here is the graph of 
dy/dt = F(y). 


Figure 2.7 
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Lecture 2: A Limited-Growth Population Model 


A glance at this graph tells us that if y> 1 or>> < -1, then / =y 2 - 1 > 0, so 
solutions must increase in this region. If —1 <y< 1, then we have y' < 0, so 
solutions decrease. And if y= 1 or y = -1, we have equilibrium points. 

With this information, we know that the slope field looks as follows. 


Figure 2.8 



And our solutions behave as follows. 


Figure 2.9 


V 
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As another example, for the ODE y* = y 4 - y 2 = y 2 (y 2 - 1), we see that there 
are 3 equilibrium points: at 0, 1, and -1. The graph of the expression 
F(y) =y 2 (y 2 ~ 1) looks as follows. 


Figure 2.10 



This graph tells us what happens between the equilibria: We have y > 0 if 
\y\ > 1, while y' < 0 if 0 < \y\ < 1. So the phase line and the solution graphs 
are as follows. 


Figure 2.11 
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Lecture 2: A Limited-Growth Population Model 
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Note that there are 3 different types of equilibrium points for this ODE. The 
equilibrium point at y = 1 is called a source since all nearby solutions move 
away from this equilibrium. The equilibrium point aty = -1 is called a sink 
since all nearby solutions move closer to this point. And the equilibrium 
point aty = 0 is called a node since it is neither a sink nor a source. 


Important Terms 


carrying capacity: In the limited population growth population model, 
this is the population for which any larger population will necessarily 
decrease, while any smaller population will necessarily increase. It is the 
ideal population. 

node: An equilibrium solution of a first-order differential equation that has 
the property that it is neither a sink nor a source. 

phase line: A pictorial representation of a particle moving along a line that 
represents the motion of a solution of an autonomous first-order differential 
equation as it varies in time. Like the slope field, the phase line shows 
whether equilibrium solutions are sinks or sources, but it does so in a 
simpler way that lacks information about how quickly solutions are 


20 















increasing or decreasing. The phase line is primarily a teaching tool to 
prepare students to make use of phase planes and higher-dimensional 
phase spaces. 

sink: An equilibrium solution of a differential equation that has the property 
that all nearby solutions tend toward this solution. 

source: An equilibrium solution of a differential equation that has the 
property that all nearby solutions tend away from this solution. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap 1.3. 
Guckenheimer and Holmes, Nonlinear Oscillations. 

Hirsch, Smale, and Devaney, Differential Equations , chap 1.2. 
Roberts, Ordinary Differential Equations , chaps. 2.1 and 2.3. 
Strogatz, Nonlinear Dynamics and Chaos , chap 2.1. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , First Order Solutions, HPG 
Solver, Phase Lines, Target Practice. 


Problems 


1. Let’s begin with some ideas from calculus. 

a. What is the integral (antiderivative) of the function y(t) = t 1 + ft 
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Lecture 2: A Limited-Growth Population Model 


b. Solve the equation;; 4 - 4 y 2 = 0. 

c. For which values of t is the function y(t) = t 3 - t positive? 

d. Sketch the graph of y(t ) = t 2 - 1 . 

e. Sketch the graph of y(t) = t{2 - t). 


2. Sketch the phase line for / = 1. 


3. Sketch the phase line for / = -y. 


4. What are the equilibria for the previous differential equation? 


5. What are the equilibrium points fory = -1? 


6. What is the phase line for the differential equation y' = y\ where n is 
some positive integer? Does the answer depend on nl 


7. What is the phase line for the differential equation;;' = /(i +j)? 


8. Find all equilibrium points for the differential equation y = sin(y) and 
then sketch the phase line. 
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9. Find an example of a differential equation that has equilibria at each 
integer value and each of these equilibria is a node. 


10. Find an example of a differential equation that has exactly 2 
equilibrium points, each of which is a node. 


Exploration 


Consider the differential equation;;' =y 2 . We know what the slope field and 
phase line look like, but can you find actual solutions? How about for 
y = [y|? Notice that the phase lines are the same for these differential 
equations, but the solutions are very different. Can you find other 
differential equations that have the same phase lines but different solutions? 
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Classification of Equilibrium Points 

Lecture 3 


I n this lecture, we see that slope field allows us to completely understand 
the behavior of solutions of first-order, autonomous differential 
equations. Recall that there are 3 different types of equilibria: sinks, 
sources, and nodes. 


Figure 3.1 



sink 


Figure 3.2 



source 
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Figure 3.3 



node 


There are actually lots of different things that can happen when an 
equilibrium point is a node. For example, the differential equation / = y 2 
has a node at y = 0. But for any nonzero value of y, solutions tend away 
from 0 if the initial value y 0 > 0, whereas solutions tend to 0 if y 0 < 0. The 
differential equation / = 0 also has an equilibrium point at 0, and this 
equilibrium point is also a node, since all solutions are equilibria. And there 
are examples where an infinite collection of equilibrium points can 
accumulate on a given equilibrium point, again giving us a node. 

We can use calculus to determine whether a given equilibrium point y 0 for 
the ODE y = F(y) is a sink, a source, or a node. Since y 0 is an equilibrium 
point, we know that the graph of F crosses the y-axis at y 0 . If the graph 
of F is increasing as it passes through y 0 , then y' is positive when y > y 0 and 
y 1 is negative when y < y 0 . This says that solutions move away from y 0 
whenever y is close to y 0 , so y 0 is a source. On the other hand, if the graph 
of F is decreasing through y 0 , then similar arguments show that y 0 is a sink. 

But we know from calculus that if F'(y 0 ) > 0, then F(y) is increasing, and 
if F'(yo) < 0, then F(y) is decreasing. So this gives us the first derivative 
test for equilibrium points: Suppose we have the differential equation 
y f = F(y), and y 0 is an equilibrium point. 
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Lecture 3: Classification of Equilibrium Points 


• If F'(y 0 ) > 0, then y 0 is a source. 

• if^'(yo) < 0, then yo is a sink. 

• If F'(y 0 ) • 0, then we get no information. 

Consider y ' = F(y) = y 2 - 1. If y = 1 or y = 0, we have equilibrium points. 
Moreover, F(y) =1-2 y. This means that F'(l) = 2, so 1 is a source. 
Similarly, F'(-l) = -2, so -1 is a sink. Therefore we know that the phase 
line is as follows. 


Figure 3.4 




And our solutions behave as follows. 


Figure 3.5 


v 
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As another example, consider}/ = y 3 — y. We can factor the right-hand side 
into / = y(y + 1 )(y - 1), so we have equilibrium points aty = 0, 1, and -1. 
Because F(y) = 3y 2 - 1, we know that 


F(0) = -1, so 0 is a sink; 

F'(l) = 2, so 1 is a source; and 
F'(-l) = 2, so -1 is also a source. 


The phase line is therefore as follows, and we again know the behavior of 
all solutions. 


Figure 3.6 


Let’s take another example: y r = y 4 - y 2 . We have equilibrium points at 0 
and at +1. Since F(y) = 4y 3 - 2y, we have F'(l) = 2 > 0; so 1 is a source. 
Because F(-\) = -2, -1 is a sink. The derivative at the 0 vanishes, so we 
get no information. But what we have about the 2 equilibria at ±1 is enough 
to tell us that 0 is a node. 
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Lecture 3: Classification of Equilibrium Points 


Figure 3.7 



1 

0 



\/ 
-1 o 
/ V 


We now turn to perhaps the most important theorem regarding first-order 
differential equations, the existence and uniqueness theorem. Suppose we 
have a differential equation / = F(y), and suppose the function F(y) is nice 
at the point y Q . (For our purposes, “nice” means that the expression is 
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continuously differentiable in y at the given point y 0 .) Then the theorem 
states that 


• there is a solution to the ODE that satisfies y(t 0 ) = yo , and 

• that this is the only solution that satisfies y(t 0 ) = y 0 . 


Let’s work a couple of examples. First, / = sin 2 (y). We have equilibrium 
points at y = 0, ± n, ± In,.... At all other points, the slope field is positive. 
So all other solutions are always increasing. However, any solution that 
starts between y = 0 and y = ;rcan never cross these 2 lines. The solution is 
trapped in this region, and since it is always increasing, the solution must 
tend to y = 0 in backward time and to y = n in forward time. We have 
similar behavior between any 2 adjacent equilibrium points. 


Figure 3.8 
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Lecture 3: Classification of Equilibrium Points 


Now let’s consider y' = ylt. Here we have problems when t = 0. The right- 
hand side of the ODE is not even defined at t = 0, never mind continuously 
differentiable there. It is easy to check that y(t) = kt is a solution of this 
ODE for any value of k. For each such k , we have y(0) = 0, so we now have 
infinitely many solutions that satisfy y(0) = 0. But there are no solutions that 
satisfy y( 0) = y 0 for any nonzero value of y 0 . Even worse, when the right- 
hand side of the ODE is not nice, the numerical methods we are using to 
plot solutions may fail. For example, the differential equation 


is not defined at y = 1. Look at what happens when we try to use the 
computer to plot the corresponding solution. 


Figure 3.9 
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Something is clearly wrong here; the numerical method the computer is 
using has failed. We will describe why this happens in Lecture 6. 
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Important Term 


first derivative test for equilibrium points: This test uses calculus to 
determine whether a given equilibrium point is a sink, source, or node. 
Basically, the sign of the derivative of the right-hand side of the differential 
equation makes this specification; if it is positive, we have a source; 
negative, a sink; and zero, we get no information. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 1.5 and 1.6. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 7.2. 

Roberts, Ordinary Differential Equations , chap. 2.2. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 2.5. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , HPG Solver, Phase Lines. 


Problems 


1. First let’s do some calculus review problems. 

a. For which values of t is the function y(t) = 1 + t 2 increasing? 

b. For which values of t is the function y(t) = \t\ decreasing? 

c. What are the roots of the function y(t) = t 2 + 2t + 1? 

d. Where is the previous function increasing? 
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Lecture 3: Classification of Equilibrium Points 


e. Sketch the graph of y(t) = t 2 + 2t + 1. 

f. Sketch the graph of y(t) = sin(f). 

g. Compute the derivative of cos(3£ + 4). 


2. Find the equilibrium points for y' =y + 1, and determine their type. 


3. Find the equilibrium points for/ = -y 2 , and determine their type. 


4. Sketch the phase line for / = -y 2 . 


5. Find the equilibrium points for / =y 2 - 1 and determine their type. 


6. Classify the types of equilibrium points for / = y 3 - 1. 


7. Does the existence and uniqueness theorem apply to the differential 
equation / = [y|? 


8. Consider / = y 2 - A, where A is a parameter. Find all equilibrium 
points and determine their type (your answer will, of course, depend 
onv4). 



9. Consider y' = Ay( 1 - y), where A is a parameter. Find all equilibrium 
points and determine their type. 


10. What is the actual behavior of solutions of 


i+F 

1 - 


Exploration 


We saw that the differential equation in problem 10 causes problems when 
we try to solve it numerically. Can you come up with other examples of 
first-order ODEs that break the numerical algorithm your computer uses to 
solve them? 
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Bifurcations—Drastic Changes in Solutions 

Lecture 4 


W e now introduce the concept of bifurcation, one of the subthemes 
in this course that will arise over and over. To bifurcate basically 
means to split apart or to go in different directions. In the theory 
of differential equations, to bifurcate means to have solutions suddenly veer 
off in different directions. In order to have such changes, our differential 
equation must depend on some sort of parameter. Sometimes, when we 
change this parameter just a little, we find a major change in the behavior of 
solutions; this is a bifurcation. 

Let’s start with a simple example: / =y 2 + A. We know that for A > 0 there 
are no equilibrium points; for A = 0 there is 1 equilibrium point, a node; and 
for A < 0 there are 2 equilibrium points, at +(— A) 1/2 . Because F'(y ) = 2 y, by 
first derivative test, the positive equilibrium point is a source and the 
negative is a sink. Look at the function graphs. 


Figure 4.1 



We see a bifurcation at ^4 = 0. 
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In our limited population growth model 
y' = ky( 1 ~y/N) 

we have 2 parameters: the growth constant k and the carrying capacity N. 
If we vary these parameters, not much happens to our solutions (as 
long and both k and N remain positive). Now let’s introduce a new 
parameter into this model to account for harvesting. Suppose our 
population is a fish population that is subject to a certain amount of fishing 
governed by the amount of fishing licenses that have been issued. 
We assume that the rate of growth of this fish population goes down 
depending on the amount of licenses issued, which is essentially the 
harvesting rate. Let h be this harvesting rate, so h > 0. Then one differential 
equation that illustrates this is the limited population growth model 
with harvesting: 

y = ky( 1 -y/N) - h. 

For simplicity, let’s revert to our easier model: 
y'=y(\-y)-h. 

First let’s find the equilibrium solutions of this equation. Look at the graph 
of the right-hand side of the ODE (i.e., the graph of y( 1 - y) - h) to see 
where this graph crosses the y-axis. When h = 0, there are only 2 places 
where y f = 0, namely y = 0 and y = 1. But as h increases, the graphs move 
downward and the 2 roots move closer together. At some h-\ alue, they 
merge into a single root. Then, for slightly higher values of h , there are no 
roots at all. So our equilibrium points have disappeared. This is the drastic 
change in the behavior of solutions—this is our bifurcation point. 
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Lecture 4: Bifurcations—Drastic Changes in Solutions 


Figure 4.2 
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Two questions arise: When does this bifurcation occur? And what are the 
ramifications in terms of the solutions of our equation (the fate of the fish 
population)? To answer the first question, we must find the h -value for 
which our equation 

y{ 1 ~y)~h = o 

has a single root. Rewriting this equation as 

- y 2 + y ~ h = 0 

allows us to solve it using the quadratic formula. We find that the roots are 
given by 


1 ± Vl — 4 h 

q 

So there is a unique root when h = 1/4. In this case, we get a single 
equilibrium point, namely y =1/2. That is, at h = 1/4, our population levels 
off at exactly half the carrying capacity. When h < 1/4, we have a pair of 
equilibria at the 2 points q+ and q-. But when h > 1/4, these equilibria 
have disappeared. 

By looking at the graphs above, we see that the equilibrium point given by 
q+ is a sink when h < 1/4, while q- is a source. Looking at the solution 
graphs, we see that all solutions that start out above the value q- tend to the 
equilibrium solution at q+. So as long as the population is high enough to 
begin with, harvesting at this rate is fine—the population survives. Even 
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when h = 1/4, all is fine as long as our initial population starts out above the 
equilibrium point. But as soon as h goes below 1/4, disaster strikes: The 
population goes extinct. 


Figure 4.3 



^ < 1/4; population survives 



h > 1/4; population dies out 
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Lecture 4: Bifurcations—Drastic Changes in Solutions 


To summarize all of this, we can record all of our phase lines in a bifurcation 
diagram. In this plot, we let the harvesting parameter h lie along the 
horizontal axis. Over each /z-value, we superimpose the corresponding phase 
line as well as all possible equilibrium points. This enables us to see the 
dramatic bifurcation that occurs when h passes through 1/4. 


Figure 4.4 



h < 1/4 h = 1/4 h> 1/4 


This type of bifurcation is called a saddle-node bifurcation. 


As another example of a bifurcation, consider the differential equation 
y' = ay + y 3 . This equation has equilibrium points at y = 0 and y = ±4~a . 
So there is only one equilibrium point when a > 0, but there are 3 equilibria 
when a < 0. We see this easily from the graphs of ay +y 3 below. 
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Figure 4.5 


dy/dt 



a < 0 


dy/dt 



dy/dt 
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< 2>0 


By the first derivative test, we see that 0 is a sink when a < 0 and a source 
when a > 0. Meanwhile, when a < 0, the other 2 equilibria are sources. We 
see this in the bifurcation diagram below. 


Figure 4.6 


V V V 



This type of bifurcation is called a pitchfork bifurcation. 
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Lecture 4: Bifurcations—Drastic Changes in Solutions 


Important Terms 


bifurcation diagram (bifurcation plane): A picture that contains all 
possible phase lines for a first-order differential equation, one for each 
possible value of the parameter on which the differential equation depends. 
The bifurcation diagram, which plots a changing parameter horizontally and 
the y value vertically, is similar to a parameter plane, except that a 
bifurcation diagram includes the dynamical behavior (the phase lines), 
while a parameter plane does not. 

limited population growth model with harvesting: This is the same as the 
limited population growth model except we now assume that a portion of 
the population is being harvested. This rate of harvesting can be either 
constant or periodic in time. 

saddle-node bifurcation: In an ODE, this is a bifurcation at which a single 
equilibrium point suddenly appears and then immediately breaks into two 
separate equilibria. In a difference equation, a fixed or periodic point 
undergoes the same change. A saddle-node bifurcation is also referred to as 
a tangent bifurcation. 

pitchfork bifurcation: In this bifurcation, varying a parameter causes a 
single equilibrium to give birth to two additional equilibrium points, while 
the equilibrium point itself changes from a source to a sink or from a sink to 
a source. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 1.7. 
Guckenheimer and Holmes, Nonlinear Oscillations , chap. 3.1. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 1.3. 
Strogatz, Nonlinear Dynamics and Chaos , chap. 3. 
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Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Phase Lines. 


Problems 


1. Find the equilibrium points of the differential equation;;' = y - A, and 
determine their type. 


2. Are there any bifurcations in the previous differential equation? 


3. Find the equilibrium points of the differential equation / = A, and 
determine their type. 


4. At which ^-values does a bifurcation occur in y' = A? 


5. Find the equilibrium points for / = y + Ay 1 . 


6. Describe the bifurcation that occurs in the family y f = Ay. 


7. Plot the bifurcation diagram for y'= Ay- y 2 . 
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8. At which type of equilibrium points would you expect a bifurcation to 
possibly occur in the equation in problem 7? 


9. Consider more generally y' = B + Ay-y 2 . Fix a value of B , and plot the 
bifurcation diagram that depends only on A. For which values of B do 
you see a different type of picture? 


10. Consider / = B + Ay - y 3 . Fix a value of B, and plot the bifurcation 
diagram that depends only on A. For which values of B do you see a 
different type of picture? 


Exploration 


Unlike in our example in this lecture, harvesting is not always constant. For 
example, harvesting of fish in certain locales is seasonal. Hence it makes 
more sense to allow the harvesting term to be a periodic function. Toward 
that end, use a computer to investigate the limited population growth model 
with periodic harvesting given by 

y =y( 1 ~y)~ h( 1 + sin(2 rt)). 

Explain the behaviors and bifurcations that you see. 
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Methods for Finding Explicit Solutions 

Lecture 5 


I n this lecture we describe several standard methods for solving certain 
types of first-order differential equations, namely, linear and separable 
first-order equations. To find these explicit solutions, we need to invoke 
another tool from calculus, integration (or antidifferentiation). Finding an 
(indefinite) integral or antiderivative is basically the opposite of finding the 
derivative; the integral of a given function is just a function whose 
derivative is equal to the given function. 

For example, the integral of the function t n for n > 0 is the function 
t n+l /(n + 1) + C, where C is some constant, since the derivative of 
t n+l /(n + 1) + C is (n + \)t n /(n + 1) = t. We denote the integral of the 
function y(t) by 

\y(t)dt. 


The dt indicates that the independent variable is t. 

Other examples from calculus that we will use are the following, in which 
ln(£) is the natural logarithm function. 

j> dt = e la /k + C 

f — = ln|/|+C 

J t 

f— = -ln|l-f|+C 
•* l-t 
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Lecture 5: Methods for finding Explicit Solutions 


A first-order linear differential equation (with constant coefficients) is a 
differential equation of the form 

y' + ky = G(t). 

The function G(t ) may be nonlinear, however, and this expression is 
sometimes called the forcing term. If G(t ) = 0 (i.e., we have no forcing 
term), then our equation is simply / = ~ky. This is called a first-order linear 
homogeneous equation, and we know its general solution, y(t) = Ce~ kt . To 
solve these equations, we use the “guess and check” method. We know how 
to solve the homogeneous part of this equation; the solution is y(t ) = Ce~ kt . 
We then make an appropriate guess to find the general solution of the 
nonhomogeneous equation. 

As an example of this method, consider y' + y = e 3t . The solution of the 
homogeneous equation y f + y = 0 is y(t) = Ce So we would make a guess 
of the formy(^) = Ce + Ae 3t . The question is, what is A1 Plugging this into 
the ODE yields A = 1/4, so our specific solution isy(f) = Ce~ l + (l/4)e 3t . The 
solution resembles e 1 as time goes backward and e 3t as time goes forward. 
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Another example of a first-order linear ODE is Newton’s law of cooling. 
Suppose you have a cup of coffee whose temperature at the moment is 
170°. Suppose also that the ambient temperature is 70°. Can you find a 
function y(t) that gives the temperature of your coffee as a function of ft 
Newton’s law of cooling says that the rate of cooling is directly 
proportional to the difference between the current temperature and the 
ambient temperature. As an ODE, this says that / = k(y - 70). If we assume 
that the coffee is initially cooling at a rate of 20° per minute, our equation at 
time t= 0 reads -20 = /(0) = k(y( 0) - 70). So k = -0.2, and the differential 
equation is / = -0.2 (y - 70) or / + 0.2 y = 14, which is first-order, linear, 
and nonhomogeneous. 

We know how to solve the homogeneous equation / = -0.2 y; the general 
solution of this equation is y(t) = Ce~ 02t . So plugging this into the left-hand 
side of the ODE yields 0. We want to find a function that we can plug into 
the left-hand side to get out the constant 14. So let’s guess a solution of the 
form y(t ) = Ce~ 02t + A, where A is some constant. Plugging this expression 
into the differential equation yields -0.2 Ce~ 02t + 0.2 Ce~ 02t + 0.2 A = 14, so 
A = 70. Therefore our explicit solution is Ce~ 02t + 70. We want the solution 
that satisfies y(0) = 170, so we put this initial value into our solution, which 
yields 170 = Ce° + 70, or C = 100. Thus our specific solution is 
y(t) = I00e~° 2t + 70. 

Another method for solving this equation is separation of variables. This 
equation is separable because we can get all the/s on the left and all the f s 
on the right. That is, we can rewrite the differential equation given by 

^ = -0.2(7-70) 

at 

in this manner: 


y-10 
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Now we integrate both sides of this equation, the left side with respect to y 
and the right with respect to t : 



The integral of the left is just In (y - 70) + constant while on the right we 
find -0.2 1 + constant. Lumping both of the constants together and calling 
them D, we find that 

ln(y-70) = -0.2f + Z>. 


Exponentiating both sides, we find 


,— 0.2 t+D ^ — 0.2 1 


y-10 = e 


where we have written the term e D as the constant C. Thus, exactly as 
above, we find the solution y(t) = Ce~ 02t + 70, and this is easily seen to be 
the general solution. 


Now let’s return to the limited population growth model given by 
y' = y( 1 ~ y). This equation is separable, so we are left with doing the 
2 integrals: 
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The right-hand integral yields t + constant, whereas the left-hand integral is 
more complicated. We can break up this complicated fraction into 2 “partial 
fractions” this way: 

__L__ 1 1 

y(l-y) 7 1-y 


So our integrals on the left reduce to 


[ — + [ ^ = ln|y|-ln|l-y| + const . 

J y J 1 — y 


Let’s assume that 0 <y < 1, so both of the terms inside the absolute values 
are positive. Therefore we are left with 

ln(y) - ln(l -y) = t + C 


for some constant C. Exponentiation plus some algebra yields 


y{ 0 = 


De l 

1 + De* ' 
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Lecture 5: Methods for finding Explicit Solutions 


These functions have graphs, below, that we have seen before. 


Figure 5.2 



If we were to assume y > 1, we would find the same solution. So is this the 
general solution? Not quite. If we try to find a value of D that solves the 
initial value problem y(0) = y 0 , we must solve the equation 


To 


= y(0) = 


£>e° 

1 + De° 


D 

i+D ' 
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Doing the algebra yields 


so we have found such a D-value as long as y 0 is not equal to 1. But we 
know the solution to the initial value problem y(0) = 1; this solution is just 
the equilibrium solution y(t) = 1. So we do in fact have the general solution 
once we tack on this solution. 


Taking the limit as t goes to infinity of 

f . De* D 

y(t) = -=- 

1 + De* e~ f +D 

yields 1, so all solutions tend to the equilibrium solution y(t) = 1, just as we 
saw earlier. 


Important Terms 


Newton’s law of cooling: This is a first-order ODE that specifies how a 
heated object cools down over time in an environment where the ambient 
temperature is constant. 

separation of variables: This is a method for finding explicit solutions of 
certain first-order differential equations, namely those for which the 
dependent variables (y) and the independent variables (t) may be separated 
from each other on different sides of the equation. 
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Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 1.2. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 1.2. 
Roberts, Ordinary Differential Equations, chap. 1.3. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 2.3. 


Problems 


1. First let’s review some integral calculus. Find the integrals of the 
following functions. 

a. y{t) = t + 5 

b. y(t) = t 2 + 2t + 1 

c. y(t) = Mr 

d. y{t) = e~‘ 

e. y{t) - 0 

2. Solve the equation ln(f) -4 = 0. 

3. Solve the equation e* - 4 = 0. 
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4. Find the general solution of y' p y = 0, and sketch the graphs of 
some solutions. 


5. Find the general solution of / - y = 1, and sketch the graphs of 
some solutions. 


6 . For the unlimited population model with y > 1, check that solutions 
assume the form discussed in the lecture. 


7. Find the general solution of the ODE y f = y 2 . 


8 . Suppose you have a cup of coffee whose temperature at the moment is 
200°, and the ambient temperature is 80°. Suppose that 1 minute later 
the temperature of the coffee is 180°. Find the function y(t) that gives 
the temperature of this cup of coffee as a function of t. 


9. Find explicit solutions of the differential equation / = 1 - y 2 satisfying 
the initial conditions y( 0 ) = 0 , y( 0 ) = 1 , andy( 0 ) = 2 . 


10. What is the general solution of the equation y' = 1 ~y 2 ? 
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Exploration 


If you are familiar with other techniques of integration, you can find 
solutions of some more complicated differential equations (though most 
such equations do not have solutions that can be found explicitly). Refresh 
your knowledge of integration techniques to find the general solutions of 
the following differential equations. 


1 . y'=l +/ 

2 . y = e‘y/(l +/) 

3 . y=y+y > 
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How Computers Solve Differential Equations 

Lecture 6 


W e have used the computer several times to plot the solutions of 
differential equations. But how does a computer produce these 
solutions? Actually, there are many different numerical algorithms 
for generating these solutions. Perhaps the most widely used method is called 
Runge-Kutta 4. For simplicity, we’ll describe a much easier to understand 
though less accurate numerical method, Euler’s method. This algorithm is 
constructed in much the same way as the more accurate methods: It uses a 
recursive procedure to approximate the solutions. 

Here is the idea behind Euler’s method. We wish to approximate the 
solution to the differential equation / = F(y, t ) starting at the initial value 
(yo, to). We first choose a step size At . This step size will usually be very 
small. Then we “step” along the slope field, moving by At units in the 

/-direction at each stage. The resulting conglomeration of little straight lines 
will be the approximation to our solution. 

More precisely, starting at the given initial point (y 0 , /o), we will recursively 
construct a sequence of points (y n , t n ) for n = 1, 2, 3, ... and join each pair of 
points (y n , t n ) and (y„+ 1 , t n+ 1 ) by a straight line. This straight line is 
constructed as follows. We start at (y 0 , t 0 ) and draw the slope line out to the 
point whose /-coordinate is t 0 + At . This is the value of t x . The 
corresponding y -value over t x is y x . Then we do the same at our point 
(yi, /i): Draw the slope line and move along it to (y 2 , t 2 ), where t 2 = t\+ At . 
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Lecture 6: How Computers Solve Differential Equations 


Figure 6.1 



We therefore have recursively t n +\ = 4 + At . So the only question is how 

do we generate the value of y n +\ knowing both y n and 4 ? The answer comes 
from algebra. We have the straight line segment passing through the known 
point (y n , tn) and having slope given by the right-hand side of the ODE— 
that is to say, the slope is F(y n , 4 ), whose value we know. Therefore our 
little slope field line has equation 

y = Mt + B, 

where M = F(y m 4 ) is the slope and B is the y-intercept. To determine the 
value of B , we know that the point (y n , 4 ) lies on our line. So we have 

y n = F(y n , 4)4 + B. 


Therefore 

B yn F(y m 4)4* 


So the equation fory„+i reads 

y n+ 1 = F (yn^n)( t n +At ) + y„-F(y n ,t n )t, 

= y n + F (y n > t n )M. 
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Thus the recursive formula to generate the values y n+x and t n+ 1 is given by 


tn +1 tn A £ 


y n + 1 =y n + F{y m t n ) At . 


Below is an image of Euler’s method applied to the differential equation 
y = y 2 ~ t with inordinately large step sizes of 0.25 and 0.125. Also 
displayed is the actual solution starting at (y 0 , to). The arrows are little 
segments of the slope field. Note how when the step size is 0.25, our 
approximate solution is way off, but when we lower the step size to 0.125, 
we get a better approximation. Lowering the step size further and further 
usually gives better and better approximations. 


Figure 6.2 
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Lecture 6: How Computers Solve Differential Equations 


Perhaps the easiest way nowadays to invoke Euler’s method is to use a 
spreadsheet. Here, for example, is a spreadsheet calculation of the solution 
to the differential equation / = y - t starting from the initial position 

y(0) = 0.3. 


Figure 6.3 
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0.9 

1 

1.1 

1.2 

1.3 

1.4 

1.5 

1.6 

1.7 

1.8 
1.9 
2 

2.1 

2.2 

23 


0.3 

0.33 

0.353 

03683 


0.37513 


0.372643 


0.3599073 


0.33539803 
0.29948783 
0.24943662 
0.18438028 
0.10281831 
0.00310014 
-0.1165899 


-0.2582488 


-0.4240737 
-0.6164811 
-0.8381292 
-1.0919421 
-1,3811363 
-1.70925 


y = y -1 

Delta t« _01_ 



'2,080175 


-2.4981925 


-2,9680117 


y(1)- 0,09720272 

Euler: 0.18438028 

Error 0.08717756 


Above we have used the step size At =0.1 and displayed the 
corresponding sequence of straight line segments. 

Unfortunately, no numerical algorithm is flawless; we can always find a 
differential equation that breaks a given numerical method. For example, 
consider the differential equation;;' = eW(y). Clearly there are equilibrium 
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points at y = 0, ± n, ± 2;r,.... So a solution that starts at y{ 0) = 0.3 should 
simply increase up to y-n as time moves on. But you can see below a 

time series for what Euler’s method yields with a step size of .08. Clearly, 
the results are bad. What you are actually seeing is chaotic behavior, 
another theme we will return to later. 

Figure 6.4 



If we change the step size just a little bit to .079, we get very different 
results for Euler’s method. This is sensitive dependence on initial 
conditions, the hallmark of the phenomenon known as chaos. 
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Lecture 6: How Computers Solve Differential Equations 


Figure 6.5 



Important Term 


Euler’s method: This is a recursive procedure to generate an 
approximation of a solution of a differential equation. In the first-order case, 
basically this method involves stepping along small pieces of the slope field 
to generate a “piecewise linear” approximation to the actual solution. 
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Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 1.4 and 
7.1-7.2. 

Hirsch, Smale, and Devaney, Differential Equations , chap. 7.5. 

Roberts, Ordinary Differential Equations , chap 2.4. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 2.8. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Euler’s method. 


Problems 


1. First let’s try some refreshers from algebra in the t-y plane. 

a. What is the slope of the straight line in the plane connecting (1,1) 
to (3, 4)? 

b. What is the slope of the straight line connecting (£ 0 , yo) to (q, y\)l 

c. What is the equation of the straight line connecting (1, 0) to (0, 1)? 

d. What is the equation of the straight line connecting (1, 0) to (3, 0)? 

e. What is the equation of the straight line connecting (1, 0) to (1, 3)? 
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Lecture 6: How Computers Solve Differential Equations 


2. What is the equation of the tangent line to the graph of y = t 2 at the 
point t= 1? 


3. Consider the differential equation / = y. Given the initial condition 
to = 0, y 0 = 1 and step size 0.1, use the formula for Euler’s method to 
compute by hand t\ andyi. 


4. Repeat the above to compute t 2 and y 2 as well as t 3 and y 3 . 


5. Sketch what you have done via Euler’s method for the approximate 
solution so far. 


6. a. For the differential equation y f = y, find the solution satisfying 

f(0)=l. 

b. Then use Euler’s method with a step size of 0.1 to approximate the 
value ofy(l). 

c. Then recompute using step sizes of 0.05 and 0.01. How do these 
compare to the approximation above with step size 0.1? 

d. Using the value of e as approximately 2.718281, how does the 
error change as we vary the step size? 

7. Use Euler’s method to approximate the solution of the differential 
equation / = (1 +/)/(1 - y) with any initial condition/0) ^ 1. What 
happens here? 
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Exploration 


A slightly better numerical approximation is given by the so-called 
improved Euler’s method. This method also uses a step size At. However, 
instead of using the slope field at the point (y m 4 ) to approximate the 
solution, this method uses a line with a slightly different slope. To get this 
slope, first use the Euler’s method approximation at the point (y„, 4 ) to find 
the endpoint of the slope line at this point whose ^-coordinate is 4 + At. 
Then take the average of the slopes at these 2 points. 

This is the slope of the straight line for the improved Euler’s method. Find 
the formula for this slope line. Now use the improved Euler’s method to 
solve the above differential equation with step sizes 0.1, 0.05, and 0.01. 
How does the improved Euler’s method compare to the ordinary 
Euler’s method? 
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Lecture 7: Systems of Equations—A Predator-Prey System 


Systems of Equations—A Predator-Prey System 

Lecture 7 


W e now begin the main part of this course, which deals with 
systems of differential equations. Our first example is the 
predator-prey system. Here we assume we have 2 species living 
in a certain environment. One species is the prey population (we’ll use 
rabbits), and the other species is the predator population (we’ll use foxes). 
Denote the prey population by R(t ) and the predator population by F(t). 

We assume first that if there are no foxes around, the rabbit population 
obeys the unlimited population growth model. If, however, foxes are 
present, then we assume that the rate of growth of the rabbit population 
decreases at a rate proportional to the number of rabbit and fox encounters, 
which we can measure by the quantity RF. So the differential equation for 
the rabbit population can be written 


dR 

dt 


= aR- bRF , 


where a and b are parameters. 

For the fox population, our assumptions are essentially the opposite; so the 
differential equation for F(t) is 

dF ^ 

— = -cF + dRF . 
dt 


Again, c and d are parameters. 

Note that both of these equations depend on R(t) and F(t). This is typical of 
systems of differential equations: We have a collection of differential 
equations involving a number of dependent variables, and each equation 
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depends on the other dependent variables. So the predator-prey system is 
given by the following equations. 


JD 

—— = aR- bRF 
dt 


- = —cF + dRF 

dt 

This is a 2-dimensional system of ODEs, since we have just 2 
dependent variables. 

More generally, a system of ODEs is a collection of n first-order differential 
equations for the missing functions y u ...,y n . Each equation depends on all 
of the dependent variables y u ..., y n as well as (possibly) t. So we get 
the following. 



y n ’t) 


dv 

-r = Fn(yi,—,yn,t) 
dt 


This is an ^-dimensional system of differential equations. Later we will 
spend a lot of time looking at linear systems of equations with constant 
coefficients. A 2-dimensional version of such a system is 

x' = ax + by 

y' — cx + dy 

where a , b, c , and d are constants. 
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Lecture 7: Systems of Equations—A Predator-Prey System 


A solution of this system of ODEs is then a pair of functions, (R(t), F(t)), 
each of which depends on the independent variable. We will think of this 
pair of functions as a curve in the plane parameterized by t. We call this 
curve a solution curve. Thus the plane becomes our phase plane, in 
analogy with the phase line viewed earlier. 

Note that we know the tangent vector to any given solution curve: This is 
the vector given by the right-hand side of our pair of differential equations 
F). So the right-hand side of our equation determines a vector field in 
the phase plane. Solution curves in the R-F plane are everywhere tangent to 
the vector field. Since these vectors can often be very large (and therefore 
cross each other in ways that make viewing the vector field difficult), we 
usually scale the vector field so that all vectors have the same length; this 
scaled field is called the direction field. 


Figure 7.1 
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The point (R 0 , F 0 ) is an equilibrium point if the right-hand sides of both 
differential equations vanish at this point. This means that we have a 
constant solution for the system given by R(t) = R 0 and F(t ) = F 0 . We 
indicate this by a dot in the phase plane. So the equilibrium points for the 
predator-prey system are given by solving the pair of algebraic 
equations simultaneously: 

R' = aR~ bRF= 0, F = ~cF+dRF= 0. 

Solving this shows that there are equilibria at (0, 0) and (c/d, alb). 

Note that when R = 0, we have R' = 0, so the rabbit population stays fixed at 
0. Also, in this case, the equation for foxes becomes F' = ~cF. This is just 
our old population model, this time with population decline rather than 
growth. So, along the F-axis, we can see the phase line corresponding to 
this first-order equation. Similarly, when there are no foxes, the differential 
equation for rabbits becomes R' = aR , the unlimited population growth 
model for rabbits. Thus we see the phase line for this equation along the 
R-axis in the phase plane. 


Figure 7.2 
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Lecture 7: Systems of Equations—A Predator-Prey System 


We may then sketch the vector field in the special case where 
a = b = c = d= 1 and use numerical methods to superimpose the plots of 
various solution curves. 


Figure 7.3 


F F.R 




On the left, we see a solution in the phase plane. Here the R -axis is the 
horizontal axis and the F-axis is the vertical axis. This solution winds 
repeatedly around the nonzero equilibrium point. On the right, we have 
superimposed the corresponding graphs of both R(t ) and F(t). Now the 
£-axis is the horizontal axis, while the F- and F-axes are vertical. 

We may modify the predator-prey system by assuming that the rabbit 
population obeys the limited growth model with carrying capacity N = 1. 
For simplicity we’ll choose all the other parameters to be equal to 1, except 
d. So our new system is 

HR 

- = R(l-R)-RF 
dt 


dF 

dt 


= -F + dRF . 
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On the R- axis, our phase line is now the limited growth phase line with 
equilibrium points at R = 0 (a source) and R = 1 (a sink). On the F-axis, we 
have (as before) the fox population tending to the equilibrium at 0. Setting 
the 2 differential equations equal to zero yields 2 equilibrium points, at the 
origin and (1, 0). But if d > 1, we get a new equilibrium point at 
(R = 1 Id, F = (d ~ 1 )/d). If 0 < d < 1, there is an equilibrium point at 
(R = l/d, F = (d - 1 )/d), but the F-coordinate here is negative, so we will 
not consider this in our population model. Thus we have a bifurcation 
when d= 1 . 

But there is more to this story. When d< 1, the phase plane appears to show 
that any solution that begins with a positive fox population now tends to the 
equilibrium point at (1, 0). That is, the fox population goes extinct. But 
when d > 1 , the solutions seem to tend to the equilibrium point 
(1 Id, (d - 1 )/d), so both populations seem to stabilize as time goes on. This 
is another example of a bifurcation. 


Figure 7.4 




Phase plane for d = 2 


Phase plane for d = 1/2 
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Lecture 7: Systems of Equations—A Predator-Prey System 


Many new types of solutions appear in systems of ODEs. For example, 
there may be periodic solutions as in the original predator-prey system. 


Figure 7.5 



And solutions may spiral toward or away from such a periodic solution. For 
example, consider the following system of equations. 

x' = (1 — yjx 1 + y 2 )x-y 
y' = x + ( 1-^/x 2 +y 2 )y 


It is very rare to be able to find periodic solutions explicitly for systems of 
ODEs, but here we can. Note that when x 2 + y 2 =1, the system reduces to 
x' - —y and y f = x. Do you know a pair of functions that satisfies this 
relation? Sure—from trigonometry we know that cos 2 (f) + sin 2 (f) = 1, and 
from calculus we know that the derivative of cos(f) is while the 

derivative of sin(f) is cos(f). So x(t) = cos(7), y(t) = sin(7) is a solution of this 
system, which is the periodic solution we see above. 
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Solutions spiral in toward a periodic solution (below). 


Figure 7.6 



Another type of differential equation that comes up often is a second-order 
equation. This is an equation that involves bothy' andy". We will deal for 
the most part with second-order equations of the form 

Y' + ay ' + by = F(t), 

where F(t ) is considered a forcing term. For example, in the next lecture we 
will consider the mass-spring system given by y" + by' + ky = 0. We can 
write this equation as a system by introducing the new variable v = y' 
(the velocity). Then our system becomes 

y = v 

v f = ~ky - bv. 

Note that this is a linear system of ODEs. This will enable us to see 
solutions both as graphs of y(t) and also as curves (y(t ), v(£)) in the 
phase plane. 
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Lecture 7: Systems of Equations—A Predator-Prey System 


Important Terms 


direction field: This is the vector field each of whose vectors is 
scaled down to be a given (small) length. We use the direction field 
instead of the vector field because the vectors in the vector field are 
often large and overlap each other, making the corresponding solutions 
difficult to visualize. Those solutions and the direction field appear within 
the phase plane. 

phase plane: A picture in the plane of a collection of solutions of a system 
of two first-order differential equations: x' = F(x , y) and y' = G(x, y). Here 
each solution is a parametrized curve, (x(£), y(t)) or (y(t\ v(0). The term 
“phase plane” is a holdover from earlier times when the state variables were 
referred to as phase variables. 

predator-prey system: This is a pair of differential equations that models 
the population growth and decline of a pair of species, one of whom is the 
predator, whose population only survives if the population of the other 
species, the prey, is sufficiently large. 

solution curve (or graph): A graphical representation of a solution to the 
differential equation. This could be a graph of a function y{t) or a 
parametrized curve in the plane of the form (x(f), y(f)). 

vector field: A collection of vectors in the plane (or higher dimensions) 
given by the right-hand side of the system of differential equations. Any 
solution curve for the system has tangent vectors that are given by the 
vector field. These tangent vectors (and, even more so, the scaled-down 
vectors of a corresponding direction field) are the higher-dimensional 
analogue of slope lines in a slope field. 
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Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 2.1. 
Edelstein-Keshet, Mathematical Models in Biology , chap. 6.2. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 11.2. 
Roberts, Ordinary Differential Equations, chap. 10.5. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Predator Prey. 


Problems 


1. Find all of the equilibrium points for the system below. 

x' = x+y 

y=y 

2. Find the equilibrium points for the system below. 

x' = x +y 
y = x+y 
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3. a. Sketch the direction field for the system below. 

x' = x 

y'=y 

b. What happens to solutions of this system? 

c. Find a formula for the solutions of this system. 

4. Write the second-order differential equation y" = y as a system of 
differential equations. 


5. a. Sketch the direction field for the following system of 

differential equations. 

x-y 

y — _ x 

b. Can you find some pairs of solutions (x(t), y{t)) of this system of 
differential equations? 

6. a. Sketch the direction field for the following system of 

differential equations. 

x' = x(l -x) 

y^ry 

b. What are the equilibrium points for this system? 

c. This system decouples in the sense that the equation for x' does not 
depend upon y, and the equation for x' does not depend on x. Use 
this fact to describe the behavior of all solutions in the phase plane. 
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Exploration 


It is impossible to find explicit solutions of the predator-prey equations 
described in this lecture. However, there is a technique from multivariable 
calculus that can give explicit formulas for where solutions lie. That is, 
we can find a real valued function H(x, y) that is constant along 
each solution. So plotting the level curves of this function shows us the 
layout of the solutions. Your job is to use calculus to find such a function. 
Hint 1: Find a function H such that dH/dt is identically equal to zero. 
This will involve the chain rule. Hint 2: A function of the form 
H(x , y) = F(x) + G(y) works. 
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Second-Order Equations—The Mass-Spring System 

Lecture 8 


I n this lecture, we introduce a slightly different kind of differential 
equation, second-order (autonomous) linear differential equations. 
These are differential equations of the form 

y" + ay' + by = G(t). 


As usual, we begin with a model, this time the mass-spring system, also 
known as the harmonic oscillator. Imagine holding a spring hanging with a 
weight (the mass) attached. If you pull the mass downward (or push it 
straight upward) and let it go, the mass will then move along a vertical line. 
Let the position of the mass be y = y(t). The place where the mass is at rest 
would be called y = 0, whereas y < 0 if the spring is stretched and T > 0 if 
the spring is compressed. 

To write the differential equation for the motion of the mass, we invoke 
Newton’s law that says that the force acting on the mass is equal to the mass 
times its acceleration. There are 2 types of forces acting on the mass: One is 
the force exerted by the spring itself, and the second is the force arising 
from friction (like air resistance). For the force exerted by the spring, we 
invoke Hooke’s law that says that the force exerted by the spring is 
proportional to the spring’s displacement from its rest position and is 
exerted toward the rest position. So this force is —ky. Here the constant k > 0 
is the spring constant. 

For the force arising from friction, we make the simple assumption that the 
force is proportional to the velocity. So this damping force is given by ~by\ 
where b is the damping constant. Here we either have b = 0 (the undamped 
case) or b > 0 (the damped case). The minus sign here indicates that the 
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damping pushes against the direction of the motion, thereby reducing the 
speed. This gives us the equation for the damped harmonic oscillator: 

my” = -by' - ky. 

To simplify matters, we usually assume that the mass equals 1. So the 
equation for the damped harmonic oscillator is 

y” + by' + ky = 0. 

Usually such a second-order differential equation comes with a pair of 
initial conditions: the initial position of the spring y( 0) and its initial 
velocity j/(0). 

We can write the equation for the harmonic oscillator as a system of 
differential equations by introducing a new variable, the velocity v = y f . So 
the system of equations is 

y = v 

V = ~ky - bv, 

where k > 0 and b > 0. Since k > 0, we see that the only equilibrium point 
for this system is at the origin (i.e., where the mass attached to the spring is 
at rest). We now have 3 different views of the mass-spring system: the 
actual motion of the mass on the left, the solution in the phase plane in the 
center, and the graphs of both position y(t) and velocity v(t) on the right. 

We also get 3 different types of solutions: one where the mass oscillates 
back and forth as it goes to rest, one where it oscillates forever, and one 
where it proceeds directly to the rest position. 
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Figure 8.1 
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How do we solve such a second-order linear equation? Recall that, for the 
analogous first-order equation y' = ky , we always had the solution 
y(t) = Ce kt . So why not guess a solution of the form e st and try to determine 
what value of s works? Plug e st into the equation y" + by' + ky = 0. We find 

e st (s 2 + bs + k) = 0. 

Therefore any root of the equation s 2 + bs + k = 0 yields an exponential 
function that solves the equation. This quadratic equation is called the 
characteristic equation; we will see many other cases where this equation 
arises. By the quadratic formula, we find that the roots are 

-Z)±VZ) 2 -Ak 
2 

But there is a problem here. What happens if the quantity b 2 - 4k is 
negative? Then the roots of our characteristic equation are complex 
numbers. What is the exponential of a complex number? We’ll tackle that 
problem a little later. Notice also that, if b 2 - 4k = 0, then we get only 1 root 
of the characteristic equation, namely —bl 2. So we get 1 exponential 
solution e ht/1 . This will not give enough solutions to solve every initial 
value problem, so more must be done here as well. We’ll come back to this 
situation again later. 

Let’s suppose that b 2 - 4k > 0. Then we have 2 distinct and real roots of the 
characteristic equation. Let’s call them a and /?, where the following 

are true. 

a = (-b-^lb r ^Ak)/2 
p = {-b + 4b r ^Ak)!2 
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Clearly, a < 0 . We also have j3 < 0 since b 2 - 4k <b 2 (remember that k , 
our spring constant, is always positive). So we have 2 decreasing 
exponential solutions, exp (a t) and exp (/? t ). And you can check that the 
general solution is 

k x e at +k 1 e pt . 

Anytime we find 2 such real, distinct, and negative roots of the 
characteristic equation, our harmonic oscillator is overdamped. This means 
that the mass slides down to rest without oscillating back and forth. 


As an example, suppose the spring constant is 2 and the damping constant is 
3. So the equation for this mass-spring system is 


y" + 3/ + 2y = 0. 


The characteristic equation is 

s 2 + 3s + 2 = 0, 

whose roots are -1 and -2. So the general solution here is 

y(t ) = k\e l + k 2 e~ 2t . 
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Here are the phase plane and the y(t) plots for several of these solutions. 


Figure 8.2 



Note that both of these solutions tend directly to the equilibrium point at the 
origin. The corresponding motion of the mass-spring system is not the usual 
oscillation back and forth, eventually coming to rest. It’s as if we were 
playing with this mass-spring system in an environment with a lot of 
resistance, like under water. This is the behavior of the overdamped 
harmonic oscillator. 

So let’s now turn to the undamped harmonic oscillator. So b = 0, and the 
equation is y" + ky = 0 or y” = —ky. Do you know a function for which the 
second derivative is —k times the function? What if k = 1? Do you know a 
function for which y" = —yl Of course, both sin(7) and cos(f) have this 
property. For other ^-values, sm(k m t) and cos {k m t) have this property. And 
we can easily check that the general solution is of the form 

c i sin (k m t) + c 2 cos (k l/2 t). 

That’s why we see the periodic behavior of solutions; any such function is 
periodic with period 2 nlk m . 
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Figure 8.3 



Important Terms 


characteristic equation: A polynomial equation (linear, quadratic, cubic, 
etc.) whose roots specify the kinds of exponential solutions (all of the form 
e At ) that arise for a given linear differential equation. For a second-order 
linear differential equation (which can be written as a two-dimensional 
linear system of differential equations), the characteristic equation is a 
quadratic equation of the form X 2 - T1 + D, where T is the trace of the 
matrix and D is the determinant of that matrix. 

damping constant: A parameter that measures the air resistance (or fluid 
resistance) that affects the behavior of the mass-spring system. Contrasts 
with the spring constant. 

mass-spring system: Hang a spring on a nail and attach a mass to it. Push or 
pull the spring and let it go. The mass-spring system is a differential equation 
whose solution specifies the motion of the mass as time goes on. The 
differential equation depends on two parameters, the spring constant and the 
damping constant. A mass-spring system is also called a harmonic oscillator. 

spring constant: A parameter in the mass-spring system that measures how 
strongly the spring pulls the mass. Contrasts with the damping constant. 
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Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 2.1, 3.6. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 6.2. 
Roberts, Ordinary Differential Equations , chap. 6.1. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 5.1. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Mass Spring. 


Problems 


1. Write the second-order differential equation for the mass-spring system 
with mass = 1, spring constant = 2, and damping constant = 3. 


2. Convert this second-order equation to a system. 


3. Find all equilibrium points for this system. 


4. Compute the second derivative of the functions sin(2£) and cos(2f). 


5. Sketch the curve ( e \ e ( ) in the plane. 
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Lecture 8: Second-Order Equations—The Mass-Spring System 


6. For the differential equation y" + 3 y' + 2y = 0, show that 
y(t) = k\e f + k 2 e~ 2t is the general solution. Do this by showing that you 
can generate a solution of this form that solves any initial value 
problem y(0) = A, /(0) = B. 


7. a. Find the general solution of the mass-spring system with spring 
constant 6 and damping constant 5. 

b. Find the solution of the above mass-spring system that starts at 
y(0) = 0 with y'(0) = 1. 

c. Sketch both the graph of y(t) and the corresponding curve (y(0, 
v(f)) in the phase plane. 

d. Explain in words what actually happens to this mass as 
time unfolds. 


Exploration 


Later in the course, we will investigate what happens to the mass-spring 
system when we apply an external periodic forcing. As a prelude to this, 
investigate what happens when we apply other, easier forcing terms. 
Consider first the mass-spring system in problem 2 above, with spring 
constant 6 and damping constant 5. What happens if we apply a constant 
force, say 1, to this system? That is, what happens to solutions of the ODE 
y" + 5 y' + 6y= 1? How about the systems where the forcing is 1 from time 0 
to time 5, and then we turn the forcing off? Find the solution to this 
equation that satisfies y(0) = /(0) = 0. 
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Damped and Undamped Harmonic Oscillators 

Lecture 9 


N ow we are ready to complete our investigations of the mass-spring 
differential equation y” + by' + ky = 0. In the last lecture, we saw 
that there is an exponential solution of the form e st when s is a root 
of the characteristic equation s 2 + bs + k = 0. But this quadratic equation 
may have complex roots, so the question is what happens to solutions when 
this occurs? Keep in mind that we also saw that the general solutions of the 
undamped equation y" +y = 0 involved sines and cosines. 

So how did we end up with trigonometric solutions to this differential 
equation instead of exponentials? Our differential equation is y" + y = 0, and 
a guess of e st yields the characteristic equation s 2 + 1 =0. The roots of this 

equation are the imaginary number i - V^l and its negative, so we have a 
solution of the form e lt . What is this complex expression? Recall from 
calculus that we can express certain functions using their power series 
expansions. That is, assuming our function f(t) is infinitely differentiable at 
t = 0, we can write/(O as below. 

fit)=m+f \o)t+/ "(o) f+/ (3) (o) f+/ <4) (o) f+• • • 


In the case of e lt , we can write this expression in a power series as below. 


it t . (it) 2 (it) 3 (it) 4 (it) 5 

2! 3! 4! 5! 

, . t 2 .t 3 t 4 J 5 

2! 3! 4! 5! 
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Lecture 9: Damped and Undamped Harmonic Oscillators 


Collecting the terms that involve the constant i and those that do not yields 


„ , t 1 t 4 t 6 ( t 3 t 5 f 

• — 1-1-h ... + l t -1-h ... 

! 4! 6! I 3! 5! 7 


But look at the right-hand side here: The first power series is just that of the 
cosine function, while the second series is the sine function. That is, we 
have one of the most famous and beautiful of all formulas in mathematics, 

Euler’s formula. 

e lt = cos(7) + i sin(7) 


Letting t = n in this formula, we get every mathematician’s favorite 
formula: e l7t +1 = 0. 

Now back to the differential equation y" + y = 0. We have the 
complex solution 

e lt = cos(f) + /sin(f). 

But our differential equation was a real differential equation; there were no 
complex numbers in sight. This means that both the real part of our 
complex solution, cos(7), and the imaginary part of the solution, sin (t), are 
also solutions. This is what we saw before: Our general solution was the 
combination of the real and the imaginary parts of our complex solution: 

y(t) = k\COs(t ) + A: 2 sin(^). 

More generally, for the differential equation y' + ky = 0, we get the 
general solution 

y(t) = k x cos (yjkt) + k 2 sin (yfkt), 

whose behavior is also periodic. 
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Figure 9.1 



Now consider the mass-spring system with spring constant k = 1.01 and 
small damping constant b = 0.2. The characteristic equation is s 2 + 0.2s + 
1.01 =0, and the roots are the complex numbers -0.1 + i and -0.1 - i. So 
we have a complex solution of the form y{t) = We can use 

properties of the exponential function and write it in the form y(f) = e~ 0At e lt . 
So by Euler’s formula, we now have complex solutions given by 
e~° At (cos(t) + zsin(f)). Taking the real and imaginary parts of this solution 
gives the general solution 

kie~ ou cos(t) + k 2 e~ 01t sm(t). 

Now we have a decaying exponential as part of our solution. So the solution 
does oscillate around its rest position, but now the amplitude of the 
oscillation (the up and down height of the y(t) graph) tends to zero as time 
goes on. We see that the mass does tend to its rest position, but it oscillates 
back and forth as it does so. This is the underdamped harmonic oscillator. 
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Lecture 9: Damped and Undamped Harmonic Oscillators 


Figure 9.2 
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More generally, when we have complex roots of the characteristic equation, 
say a + ib and a - ib (here a is necessarily negative), we get the 
general solution 

k x e at cos (bt) + k 2 e at sin(Z^), 
which similarly oscillates down to its rest position. 


Recall that there was one other case of the harmonic oscillator that we have 
not yet solved, the case for which we found only one root for the 
characteristic equation. This is the critically damped mass-spring case. 
For example, consider the mass-spring system 

y" + 2/ +y = 0. 
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The characteristic equation is now s 2 + 2s + 1, which factors into 
( s + l) 2 = 0. So we have a pair of real roots given by s = —1. We know one 
solution, namely, e~\ As before, we can multiply this expression by any 
constant to find other solutions of the form k { e \ What about other 
solutions? We can easily check that y(t ) = te~ f also solves this differential 
equation and that, moreover, 

kie~ f + k 2 te t 

is the general solution. What happens to the solution given by te~ f as time 
goes on? Using PHopital’s rule, we find that the limit of the expression 
te~ f also tends to zero. Thus we see that both terms in the general solution 
go to zero—and in fact, in the phase plane, they do so just as in the 
overdamped case. 


Figure 9.3 



Vo = 


v 0 = 


d 2 y/dt 2 = -(k/m)y - (b/m)v 
dy/dt = v 

dv/dt = -(k/m)y - (b/m)v 
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Lecture 9: Damped and Undamped Harmonic Oscillators 


Important Terms 


Euler’s formula: This incredible formula provides an interesting 
connection between exponential and trigonometric functions, namely: The 
exponential of the imaginary number (i-t) is just the sum of the 
trigonometric functions cos(7) and i sin(7). So e^ li) = cos (t) + i sin(7). 

critically damped mass-spring: A mass-spring system for which the 
damping force is just at the point where the mass returns to its rest position 
without oscillation. However, any less of a damping force will allow the 
mass to oscillate. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 2.1 and 3.6. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 3.4. 

Roberts, Ordinary Differential Equations , chap. 6.1. 

Strogatz, Nonlinear Dynamics and Chaos , chaps. 5.1 and 7.6. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Mass Spring, RLC Circuits. 


Problems 


1. Use the product rule to compute the derivative of e\os(3 1). 


2. Use Euler’s formula to expand e^ 2t + 3lt \ 
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3. Solve the quadratic equation s 2 + bs + k = 0, where b and k 
are constants. 


4. The point (sin(7), cos(7)) lies on the unit circle in the plane. In which 
direction does it move as t increases? 


5. Is the mass-spring system y" + 5/ + 6y overdamped or underdamped? 


6. Use 1’HopitaTs rule from calculus to verify that the limit as t 
approaches infinity of te~* is 0. 


7. Find the general solution of the mass-spring system with spring 
constant 2 and damping constant 2. 


8. If we start the mass-spring system in problem 7 off from its rest 
position but with velocity equal to 1, what is the maximum distance the 
center of mass will move? 


9. Consider the mass-spring systems with spring constant 1 and damping 
constant b. For which values of b is this system underdamped, 
overdamped, undamped, or critically damped? 
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Lecture 9: Damped and Undamped Harmonic Oscillators 


10. We know the general solution of the undamped harmonic oscillator 
given by y" +y = 0. What happens if we apply a constant force equal to 
1 to this system? That is, what is the general solution of the equation 
y" + y = 1? And what is the motion of the corresponding mass-spring? 


Exploration 


One of the simplest examples for electric circuit theory is the RLC circuit, 


Figure 9.4 


R 

|—WV 

W© 


L 




where R is the resistance, L is the inductance, C is the capacitance, and V^t) 
denotes the voltage at the source. The differential equation for the voltage 
across the capacitor v(t) is similar to our mass-spring example and is 
given by 


LCv" + RCv' + v=Vs(t). 

Find the solution to this equation when V s (t) is zero. Compare the behavior 
of these solutions to those of the mass-spring system. 
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Beating Modes and Resonance of Oscillators 

Lecture 10 


I n this lecture, we consider the periodically forced harmonic oscillator. 
That is, instead of hanging our spring to the ceiling, we now move the 
spring up and down periodically. So our differential equation becomes 

y" + by ' + ky = G(t\ 

where G(t ) is some periodic function. For simplicity, let us take 
G(t) = cos (wt). So the period of the forcing term is 2 tt/w. There are 2 very 
different cases to consider. The first is when the damping constant is 
nonzero. In this case, if there were no periodic forcing, our mass-spring 
system would settle down to rest. But with the periodic forcing present, the 
mass cannot end up in its rest position; rather, the mass eventually begins to 
oscillate with period 2 tt/w. 


As an example, consider the differential equation 
y" + 3/ + 2 y = cos(t). 


Solving the homogeneous equation (the equation with 0 replacing cos(£)) 
yields the solution 

y(t) = k x e~ 2 ‘ +k 2 e~‘ 

since the characteristic equation is s 2 + 3s + 2 = (5 + 2)(s + 1). So the 
unforced oscillator tends to rest just as we expected. 
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Lecture 10: Beating Modes and Resonance of Oscillators 


To find one solution of the forced mass-spring system, we need to make a 
guess. Clearly, we cannot guess just ^4sin(7), since the / term will give us a 
cosine term. Similarly, we cannot guess Bcos(t). So the appropriate guess is 
^sin(0 + 5cos(7). Plugging this guess into the equation and doing a little 
algebra yields the particular solution with A = 3/10 and B = 1/10. So any 
solution of this equation is of the form 

y(t) -k x e + + -T sin(Y) + -jLcos(Y). 

No matter which values of k\ and k 2 we take, the corresponding solution 
ends up looking like 

7(0 = ^sin(0 + -jVcos(0. 

This equation has period 2n, which is exactly the period of our forcing 
term. This is called the steady-state solution. 


Here are the graphs of 3 different solutions, together with the solution 
curves in the phase plane. Note how all solutions quickly approach the 
steady-state solution. 


Figure 10.1 
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Now let’s consider the case where there is no damping term and our forcing 
term is cos(w£), so we have a parameter present. For simplicity, consider 
the equation 

y” + y = cos(w£). 


Since we have no damping, the homogeneous equation is y" + y = 0. So the 
general solution is &iCOs(7) + & 2 sin(f). This solution does not tend to rest. To 
get the full solution, we make the guess Asm(wt) + i?cos(wf). Plugging this 
into the differential equation, we find 

B{ 1 - w 2 )cos(utf) +A( 1 - w 2 )sin(wf) =y" +y = cos (wt). 


So we must have B = 1/(1 - w 2 ) and A = 0. That is, our particular solution is 

—^—r cos (wt) + k x cos(^) + k 2 sin(^). 

1 -w 


But wait a moment—this is fine as long as w is not equal to 1; if w =1, this 
is no longer the solution. That is, if the period of the forcing is the same as 
the natural period of the spring, we are in trouble. We’ll deal with this case 
later. But first look what happens for different values of w. Here are some 
graphs of the solution that satisfies y(0) =/(0) = 0. First, when w= 1.5, we 
see a periodic motion, both in they-v plane and as the graph of y(t). 
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Lecture 10: Beating Modes and Resonance of Oscillators 


Figure 10.2 
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But when w= V2 , the solution is no longer periodic. 


Figure 10.3 
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Lecture 10: Beating Modes and Resonance of Oscillators 


And when w gets very close to 1, say w = 1.06, we get beating modes. 


Figure 10.4 
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Note that the ^-coordinate runs here from -20 to 20. So the mass sometimes 
undergoes large oscillations and sometimes does not. Why? Since w is very 
close to 1, the natural period and the forcing period are close together. This 
means that for a long period we are forcing the spring exactly the way it 
naturally wants to go, so the oscillations become larger and larger. But 
eventually, cos(t) and cos(utf) are relatively far apart, and this means that we 
are forcing in the direction opposite to that the spring wants to go, resulting 
in smaller and smaller oscillations. 


Now let’s go back to the case w = 1, that is, the differential equation 
y" + y = cos (t). 


We no longer have our old solution 

—^-cos(w0 + k x cos(f) + k 2 sin(f) 

1 - w 

since w is now equal to 1. Indeed, we cannot have a solution of the form 
Acos(t) + Bsin(t) since that expression solves the homogeneous equation. So 
what to do? How about making the guess y(t) = Atcos(t) + Btsin(t)? 
Plugging this into the differential equation then yields 

y" + y = -2Hsin(/) + 2i?cos(/). 


Setting this equal to cos (t) tells us that B = 1/2 while A = 0. So our general 
solution now becomes 

y t sin(Y) + k x cos(7) + k 2 sin(7). 
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Lecture 10: Beating Modes and Resonance of Oscillators 


Therefore the solution that satisfies y(0) =/(0) = 0 is just 




Note that this expression then oscillates back and forth with increasing 
amplitude; that is, our spring keeps moving alternately higher and lower. 
Note that in the graph below, the y-values extend from -100 to 100. 


Figure 10.5 






























This is the phenomenon called resonance. We keep pushing and pulling our 
mass-spring system in just the way it wants to go, so the oscillations keep 
getting larger and larger, eventually leading to disaster. 

Some experts blame resonance for the collapse of the Tacoma Bridge (a.k.a. 
“Galloping Gertie”) back in 1940, just 4 months after the bridge was opened 
to traffic. A similar phenomenon led to the closing of London’s Millennium 
Bridge in June 2000, just 2 days after it was opened to pedestrian traffic. 
Indeed, when crossing a bridge, soldiers no longer march in step for fear 
that their steps will coincide with the resonant frequency of the bridge. 


Important Terms 


beating modes: The type of solutions of periodically forced and undamped 
mass-spring systems that periodically have small oscillations followed by 
large oscillations. 

resonance: The kind of solutions of periodically forced and undamped 
mass-spring systems that have larger and larger amplitudes as time goes on. 
This occurs when the natural frequency of the system is the same as the 
forcing frequency. 

steady-state solution: A periodic solution to which all solutions of a 
periodically forced and damped mass-spring system tend. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 4.3 and 4.5. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 6.2. 

Roberts, Ordinary Differential Equations , chap. 6.1. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 7.6. 
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Lecture 10: Beating Modes and Resonance of Oscillators 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Beats and Resonance, 
Mass Spring. 


Problems 


1. a. Solve the periodically forced first-order equation 
y +y = sin(f). 

b. What happens to all solutions of this equation? 

c. Sketch the graphs of some of these solutions in the region 
0<t<20. 

d. What happens to solutions if we change this equation to 

y ~ y = sin(/)? 

e. What happens to these solutions as time moves backward? 


2. Find the steady-state solution of the forced mass-spring system with 
spring constant 2, damping constant 1, and forcing term sin(f). 


3. What is the general solution of the forced mass-spring system given by 

y"+y = e~ t 2 


100 







4. Consider the mass-spring system;;" + 3 y = sin(wf). For which value of 
w is this system in resonance? 


5. For which values of w is the function cos (wt) + cos (t) a 
periodic function? 


6. Consider the undamped and forced mass-spring system y" + 
y = cos(wf). Using the results of the previous problem, determine the 
values of w for which the solutions of this system are 
periodic functions. 


Exploration 


Consider the case of a pair of undamped harmonic oscillators. 

y r = ~wi 2 yi 
y 2 " = ~Wy 2 


These equations are decoupled, so we can easily solve both. But let’s view 
these solutions a little differently. The quantity w 2 /w i is called the frequency 
ratio. First show that the pair of solutions of these 2 equations is periodic if 
and only if the frequency ratio is a rational number. Plot a solution ( yi(t ), 
y 2 (0) i n the plane when, say, uq = 2 and w 2 = 5. What if w\= 1 and w 2 = 2 1/2 ? 
Change these equations to polar coordinates (r 7 , 0 7 ) for j = 1,2. What is the 
new system of differential equations? Do you see that r- = 0? So we can 
now think of solutions as residing on a torus (i.e., the surface of a 
doughnut). What would these solutions look like on the torus? 
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Lecture 11: Linear Systems of Differential Equations 


Linear Systems of Differential Equations 

Lecture 11 


W e now move on to linear systems of differential equations. 

Much of the theory surrounding linear systems of differential 
equations involves tools from the area of mathematics known as 
linear algebra, including matrix arithmetic, determinants, eigenvalues, 
and eigenvectors. 

Let’s consider 2-dimensional, autonomous linear systems. These are 
systems of differential equations of the form 

x' = ax + by 

y r = cx + dy , 

where a , b , c, and d are all parameters. We will often write this in matrix 
form Y = AY, where Y is the 2-dimensional vector 



\y) 

and A is the 2 x 2 matrix 
(a b\ 


The product of the matrix A and the vector Y is the new vector whose 
first entry is the dot product of the first row of A with Y and whose second 
entry is the dot product of the second row of A with the vector Y —that is, 
the vector 


AY = 


r ax + by^ 
K cx + dy j 
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A solution of such a linear system is then a vector depending on t (a curve 
in the plane) given by 


m= 


'm' 

^( 0 , 


where x(7) and y(t ) are the solutions of the given system. So x'(0 = ax(t) + 
by(t ), and/(/) = cx (0 + <X(0- 


For example, consider the linear system 



n 

Y = 

(i 

n 

f x l 


-b 


v3 

-b 



We can easily check that one solution is x(t) = e 2t , y(t) = e 2t . We write this 
solution in vector form as 


m- 


fe 2 '] 

2 1 

rn 


= £ 



Note that this is a straight line solution of the system; in the phase plane, 
this solution lies along the straight line y = x. It tends to the origin as time 
goes to -oo and to oo as time goes to oo . 

Another solution is x(t) = e ~ 2t , y(t) = 3e ~ 2t , or in vector form, 


e " 2 ' 1 

—2t 

r i ^ 


- e 

v-3v 
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Lecture 11: Linear Systems of Differential Equations 


Again, this is a straight line solution, this time lying on the line y =* —3x in 
the phase plane and tending to the origin as time increases. These solutions 
in the phase plane look as follows. 


Figure 11.1 
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As earlier, any expression of the form 


Y(t) = c,e 2 


(l) 

+ c 2 e 

rn 

4, 




for any values of c\ and c 2 is also a solution. Here is a collection of such 
solutions in the phase plane. 
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Figure 11.2 


V 



We say that the equilibrium point at the origin in this phase plane is a 

saddle point. 


In fact, 


Y(t) = qe 2 


fi] 

+ c 2 e 21 

rn 



v-3. 


is the general solution to this differential equation. In order to see why this 
is true, we must be able to determine values of c\ and c 2 that solve any 
initial value problem of the form 


7(0) = 


{ v \ 


U). 
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Lecture 11: Linear Systems of Differential Equations 


That is, we must be able to find C\ and c 2 so that 


Y (0) = Cj 


m m (x,' 


vly 


+ Cn. 


V-3y 


0 

uv 


Notice that the 2 vectors (1,1) and (1, -3) are linearly independent. That is, 
they do not lie along the same straight line passing through the origin. This 
means that by the parallelogram rule, we can always stretch or contract 
these vectors (by multiplying by c\ and c 2 ) so that the resulting vectors add 
up to (x 0 , To). 


Figure 11.3 
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Whenever we have 2 such solutions whose initial vectors are linearly 
independent, a similar combination of them gives us the general solution. 
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As another example, consider 


Y' = 


f-l 

vl 



Y. 


We can easily confirm that the only equilibrium solution is at the origin. To 
find the other solutions, note that the first equation reads x' = —x. So we 
know that the general solution here is k { e 1 . Then the second equation is y' = 
—2y + k x e\ which is a first-order linear and nonhomogeneous equation. We 
know how to solve this. The solution is y(t) = k 2 e~ 2t + k\e~ l . So our solutions 
to the linear system are 

x(t) = k\e l 

y(t) = k 2 e~ 2t + k\ e l . 


In vector form, this solution is as below. 


Y(t) = k x e l 


fW 


vly 


+ k,e 2t 


f 

V 


0 " 

K 


When k\ = 0 and k 2 = 1, we get the straight line solution x(t) = 0 and 
y(t) = e~ 2t lying along the x-axis. And when k\ = 1 and k 2 = 0, we get another 
straight line solution x(t) = e~\ y(t) = this time lying along the line y = x. 
It is the same as before, but now both straight line solutions tend inward to 
the origin. We see that the phase plane is a sink. 
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Figure 11.4 



Thus we see again 2 straight line solutions for this system. This will be the 
way we proceed to solve linear systems; we’ll first hunt down some straight 
line solutions and then use them to put together the general solution. 

This straight line method does not always produce straight line solutions, 
however. There are many other possible phase planes for linear systems, 
which we produce with the same method. 

For example, for the system 


we have x'(t ) = ~y(t) while y’(t) = x(t). So some solutions are x(t) = ccos(f), 
y(t) = csin(7) for any constant c. These solutions all lie along circles in the 
phase plane, and the equilibrium point at the origin is now called a center. 
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Figure 11.5 
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Let’s look at a different example. 


Y' = 


(- 0.1 


i ^ 
- 0 - 1 , 


Y 


We easily determine that some solutions are of the form 


Y(t ) = ce~ ou 


/ cos (?) ' 
v -sin(0 y 


Without the exponential term, these solutions would all lie along circles as 
above. But, since we have the exponential of-0.k, these solutions must 
now spiral in toward the origin as time increases. Our phase plane now 
looks as follows, and this equilibrium solution is called a spiral sink. 
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Lecture 11: Linear Systems of Differential Equations 


Figure 11.6 
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The first question concerning linear systems is what are the equilibrium 
points. Certainly the origin, (0, 0), is always an equilibrium point, but are 
there others? If so, they would be solutions to the system of linear 
algebraic equations 

ax + by = 0 

cx + dy = 0. 


The first equation says that x = ~(b/a)y (as long as a is nonzero). 
Substituting this into the second equation, we find 

(<ad - bc)y - 0. 
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So there are only 2 possibilities here, either y = 0 or ad - be = 0. In the first 
case, we then have that x is also 0, so we get the equilibrium point at the 
origin we already know about. In the other case, ad - be = 0, we then have 
that any y value works; so x = ~(b/a)y. So this yields a straight line of 
equilibrium solutions given by x = —(b/a)y. 

A similar argument works if a = 0, at least as long as one of the coefficients 
in the matrix is nonzero. If all the entries in the matrix equal 0, we have a 
pretty simple system—namely, x' = 0 and y' = 0, so all points (x, y) are 
equilibrium points. We’ll neglect this simple case in the future. 

So, to summarize, for the linear system of ODEs Y = AY, we have 

1. A straight line of equilibrium points if ad - be = 0 

2. A single equilibrium point at the origin if ad - be ^ 0. 


The quantity ad - be will reappear often in the sequel; it is called the 
determinant of the matrix. We abbreviate it det4 = ad - be. 


Notice also that a linear system of algebraic equations 

ax + by = 0 
cx + dy = 0 

has nonzero solutions if and only if det4 = 0. 
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Lecture 11: Linear Systems of Differential Equations 


Important Terms 


center: An equilibrium point for a system of differential equations for 
which all nearby solutions are periodic. 

determinant: The determinant of a 2-by-2 matrix is given by ad - be, 
where a and d are the terms on the diagonal of the matrix, while b and c are 
the off-diagonal terms. The determinant of a matrix^ (det4, for short) tells 
us when the product of matrix A with vector 7 to give zero (A 7=0) has 
only one solution (namely the 0 vector, which occurs when the determinant 
is non-zero) or infinitely many solutions (which occurs when the 
determinant equals 0). 

saddle point: An equilibrium point that features one curve of solutions 
moving away from it and one other curve of solutions tending toward it. 

spiral sink: An equilibrium solution of a system of differential equations 
for which all nearby solutions spiral in toward it. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 3.1. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 2.4. 
Kolman and Hill, Elementary Linear Algebra, chap. 3. 

Lay, Linear Algebra, chap. 3. 

Roberts, Ordinary Differential Equations, chap. 7. 

Strang, Linear Algebra and Its Applications, chap. 4. 

Strogatz, Nonlinear Dynamics and Chaos, chap. 5.1. 
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Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Linear Phase Portraits. 


Problems 


What is the product of the matrix 


1 3 

2 1 


and the vector 


f 




2. a. Write the system of equations 
x' = y 

y = ~ x 

in matrix form. 

b. What are the equilibrium points for the above system? 


3. Are the vectors (1,2) and (-4, -8) linearly independent? 


4. a. Is 


( „ 2t.t\ 

e+e 


a solution of the system 

x' = x +y 


y' = 2y? 
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Lecture 11: Linear Systems of Differential Equations 


b. Is 


a solution of the above equation? 


f 

e 


2 


V 


0 


/ 


5. Consider the linear system below. 

x f = x 
y' = x + 2y 

a. Write this system in vector form. 

b. Since the system can be decoupled, first find the solution of x f = x 
and then use that to find the solution of y' = x + 2 y. 

c. Write this solution in vector form. 

d. Do the 2 above solutions generate the general solution? 

e. Sketch the above solutions in the phase plane. 


Exploration 


We have seen that the determinant of a 2 x 2 matrix allows us to decide 
when a system of 2 linear equations has a nonzero solution. What about the 
case of 3 linear equations, such as the below? 

a n x+ a\ 2 y + a^z = 0 
d2\x + a2?y + CL23Z = 0 
a 3 \X + a 3 2y + CI33Z = 0 

Can you find a formula (depending on the a^) that tells you when this 
system of equations has nonzero solutions? 
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An Excursion into Linear Algebra 

Lecture 12 


Recall the linear system that we solved last time: 


Y' = 


1 1 


Y. 


There were 2 straight line solutions given by 


Y l (t) = e 2 ‘ 


f\\ 


\b 


and Y 2 {t)-e 


f i A 


v-3y 


Yi(t) lies along the diagonal line y = x, and Y 2 (t) lies along the liney = — 3x. 
Note that 



y 


f 2 ^ 

= 2 


and 

(i n 

1 1 


b 2 ^ 

= -2 

rn 

b -b 

<b 


A 


lb 


b -b 

v-sj 


U; 


I- 3 J 


We’ll see these equations again. 


How do we solve the general linear system? As above, we have often seen 
that there are solutions of a system that run along a straight line. How do we 
find these solutions? Well, we would need the vector field to either point 
straight in toward the origin or point directly away from the origin. That is, 
we would need the right-hand side of our system, namely A Y, to be some 
multiple of Y, say AY . That is, we would need to find a vector Y such 
that AY = AY or AY-AY = 0, where 0 is the zero vector—that is, (0, 0). 
Note that this is precisely what is happening in our example above. 
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Lecture 12: An Excursion into Linear Algebra 


We usually write this latter equation as {A-XI)Y - 0, where / is the 
identity matrix 


1 = 


f 1 

v0 


o' 


Written out, these equations become 

(a - X)x + by = 0 
cx + (d - X)y = 0. 


So to find places where the vector field points toward or away from the 
origin, we must find a value of X for which the above equations have a 
nonzero solution. Such a X -value is called an eigenvalue for our system, 
and any nonzero solution of the equation (A-XI)Y = 0 is called an 
eigenvector corresponding to X . 

So how do we determine the eigenvalues? We need to find a value of X for 
which the above system of algebraic equations has nonzero solutions. But 
we saw earlier that that happens when the determinant of the corresponding 
matrix is zero. So X 0 is an eigenvalue if it is a root of the equation 

det(^4 - XI) - 0. Let’s write that out explicitly. We find that 

det(^4 - XI) - (a- X)(d - X) - be = X 2 - {a + d)X + ad - be = 0. 


This equation is also called the characteristic equation. Note that it is a 
quadratic equation, so there are usually 2 distinct roots for this equation. 
These roots are the eigenvalues for the matrix. Then we can find the 
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corresponding eigenvector by simply finding a nonzero vector solution to 
the equation 


(A-A 0 I)Y = 0, 


where A 0 is the given eigenvalue. 


By the way, note that in the characteristic equation 
X 1 - {a + d)/1 + ad -be = 0, 

we see the term ad - be. That’s what we called the determinant of the 
matrix A, or det4. There is also the term a + d (i.e., the sum of the main 
diagonal terms in the matrix). This expression is called the trace of matrix 
A, which we write as Tv A. So the characteristic equation can be written 

A 2 -(Tr A)A + detA = 0. 


The important observation is that if A is an eigenvalue for the matrix A and 
7 0 is its corresponding eigenvector, then we have a solution to the system 
given by 

m = e l %, 

which is easily checked. 


For example, consider the linear system 


7 ' = 


r 2 

v0 


n 


7. 
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Lecture 12: An Excursion into Linear Algebra 


We first find the eigenvalues: They are the roots of the characteristic equation 


det 


(2-A 


1 ^ 
-l-A, 


= 0 , 


which simplifies to (2-A)(-l-A) = 0. Clearly these roots are 2 and -1. So 

these are our eigenvalues. Next we find the eigenvector corresponding to the 
eigenvalue 2. To do this, we must solve the system of algebraic equations 

(2-2)x + y = 0 
0x + (-l-2)y = 0 

that reduces to y = 0 and -3y = 0. Notice that these 2 equations are 
redundant; this is always the case since we know that we must find a 
nonzero solution to the eigenvector equation. Our solution is any (nonzero) 
vector with y-coordinate equal to 0. So for example, one eigenvector 
corresponding to the eigenvalue 2 is the vector 


(we could have chosen any nonzero x-component). 


So we have one solution to this system, namely 


m = e 2t 


a 


or x(t) = e 2t , y(t) = 0. Note that this is a straight line solution lying on the 
x-axis. It is easy to check that any constant times this solution is also a 
solution to our system, so we find other straight line solutions lying on the 
x-axis—some on the positive side and some on the negative side. But all tend 
to (plus or minus) infinity as time goes on and to 0 as time goes backward. 
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For the eigenvalue -1, we similarly find a corresponding eigenvector 


Y = 



and thereby get the solution 



This is also a straight line solution lying along the line y = ~3x, which tends 
to the origin as time goes on and away from the origin (toward infinity) as 
time goes backward. 

So we get a whole family of solutions of the form 


7 27 

k x e 

rn 


\ 1 1 

+ k 2 e 


vV 


-3, 


Indeed, this turns out to be the general solution. Below are the straight line 
solutions together with some other solutions in the phase plane. This is an 
example of a phase plane that has a saddle equilibrium point at the origin. 


Figure 12.1 


y 


v 


o- 


- 4 - 


4 - 


-4 



4 - 
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Lecture 12: An Excursion into Linear Algebra 


There are 2 other possible things that can happen when we have 2 real, 
distinct, and nonzero eigenvalues. Either both can be positive or both can be 
negative. In the former case, our equilibrium point is a (real) source, since 
our general solution will be of the form 

k/' t Y x +k 1 e^ t Y 1 , 

where both \ and are positive. Thus all nonzero solutions will move 
away from the origin. In the other case, we have a (real) sink since all 
solutions will now tend to the origin. 

For example, if our matrix is 

r 0 P 

v-1 -3/ 

then the eigenvalues are the roots of P +3/1+ 1 = 0, which by the quadratic 
-3± V5 

formula are -—-. Both of these are negative, so we get a real sink at 

the origin. 

Figure 12.2 


y 
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But if our matrix is 


r 


V 


0 

-1 


n 

3 , 


the eigenvalues are 


3±V5 

2 


, which are now both positive, so we get a 


real source. 


Figure 12.3 


V 



Then there are a number of other possibilities for the eigenvalues. These 
roots of the characteristic equation could be complex; real, nonzero, but not 
distinct (i.e., repeated roots); or zero. As we shall see in upcoming lectures, 
each of these cases gives very different phase planes. 
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Lecture 12: An Excursion into Linear Algebra 


Important Terms 


eigenvalue: A real or complex number usually represented by X (lambda) 
for which a vector Y times matrix A yields non-zero vector XY. In general, 
an n x n matrix will have n eigenvalues. Such values (which have also been 
called proper values and characteristic values) are roots of the 
corresponding characteristic equation. In the special case of a triangular 
matrix, the eigenvalues can be read directly from the diagonal, but for other 
matrices, eigenvalues are computed by subtracting X from values on the 
diagonal, setting the determinant of that resulting matrix equal to zero, and 
solving that equation. 

eigenvector: Given a matrix A, an eigenvector is a non-zero vector that, 
when multiplied by A, yields a single number X (lambda) times that vector: 
AY= XY. The number X is the corresponding eigenvalue. So, when X is real, 
AY scales the vector Y by a factor of lambda so that AY stretches or 
contracts vector Y (or —Y if lambda is negative) without departing from the 
line that contains vector Y. But X may also be complex, and in that case, the 
eigenvectors may also be complex vectors. 

trace: The sum of the diagonal terms of a matrix from upper left to 
lower right. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 3.2-3.3. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 3.1. 
Kolman and Hill, Elementary Linear Algebra, chap. 7. 

Lay , Linear Algebra, chap. 5. 
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Roberts, Ordinary Differential Equations , chap. 7. 


Strang, Linear Algebra and Its Applications , chap. 5. 
Strogatz, Nonlinear Dynamics and Chaos , chaps. 5.1-5.2. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Linear Phase Portraits. 


Problems 


1 . 


Compute the determinant of the matrix 


3^ 

i. 


2. What do you know about the equilibrium points for the system 


Y' = 


^1 

2 


3^1 


7? 


3. Compute the trace and the determinant of the matrix 

r 2 X s 

vl 2/ 


123 












Lecture 12: An Excursion into Linear Algebra 


4. a. Sketch the direction field for the linear system 


Y' = 


vO 


O''! 


Y. 


b. Do you see any straight line solutions for this system? 


5. What are the eigenvalues of the matrix 


6. What are the eigenvalues of an upper triangular matrix 


b) 

where * can be any number? 


7. What are the eigenvalues and eigenvectors for the very special matrix 

fo (A 

? 

to oj • 


8. Find the eigenvalues and eigenvectors of the matrix 

' 1 3 ^ 

v V2 3V2 y ' 
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9. Find the general solution of 


Y' = 


(3 


n 

3 , 


Y 


and sketch the phase plane. 


10. Find the solution to the previous system satisfying the initial condition 

m = (i,o). 


Explorations 


1. Consider the linear system 


Y' = 


f-l 



Y. 


First find both straight line solutions. Then determine how all other 
solutions tend to the origin as time goes on. 

2. In the last exploration of linear algebra, we saw that that the 
determinant of a 3 x 3 matrix allows us to decide when a system of 3 
linear equations has a nonzero solution. This should allow us to expand 
the notion of eigenvalue and eigenvector to this case. So find the 
eigenvalues and eigenvectors of the matrix 

r 3 0 P 
0 2 0 . 


What would be the behavior of the solutions of the corresponding linear 
system of differential equations? 
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Lecture 13: Visualizing Complex and Zero Eigenvalues 


Visualizing Complex and Zero Eigenvalues 

Lecture 13 


Recall that to solve the linear system of differential equations given by 


Y' = 


a b 


\C dj 


we followed a 4-step process. First we wrote down the characteristic equation 

X — (ci + <?)/1 + ( ad — be ) — X — T A, + D — 0 , 

where T is the trace of the matrix and D is the determinant. Second, we 
solved this quadratic equation to find the eigenvalues 


T± 


YF 2 


-AD 


K = 


Third, for each eigenvalue A ± , we found the corresponding eigenvectors V ± 
by solving the system of algebraic equations below. 

(a - X)x + by = 0 
cx + (d - X)y - 0 

Finally, this gave us the general solution below. 


k 1 e Aj V_+k 2 e A *‘V + 


But the roots of the characteristic equation may be complex—so what 
happens in that case? 
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Let’s begin with the simple example 


7 ' = 


^0 


n 

o. 


Y. 


We actually know how to solve this system since the system is just x' = y 
and y = ~x. So we really have x” = —x, and we know the solution of this 
second-order equation. The general solution is x(t) = k\COs(t) + k 2 sin(t). 
Therefore y{t) = x'(t) = -&isin(7) + & 2 cos(0. 


Altogether, this yields 


Y(t) = k i 


' cos(/) ' 
v -sin (t), 


"sin(0' 

v cos (t), 


Let’s see how this solution arises out of the eigenvalue/eigenvector process. 
Our characteristic equation here is X 1 +1 = 0, so the eigenvalues are the 
imaginary numbers i and ~i. To find the eigenvector corresponding to /, we 
must solve the system of equations 

—ix + y = 0 
—x~iy = 0. 

These 2 equations are redundant, since multiplying the second equation by i 
gives the first equation. So the eigenvector corresponding to the eigenvalue i 
is any nonzero vector that satisfies y = ix. So for example, one eigenvector is 

rn 

vL ’ 


We therefore get the solution 

■ m 

7(0 = e* 

VJ 
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Lecture 13: Visualizing Complex and Zero Eigenvalues 


But there is a problem here. We started with a system of real differential 
equations but ended up with a complex solution. Let’s remedy this using 
Euler’s formula. Our solution can be expanded via Euler to 


Y(t) = (cos(f) + i sin (7)) 




r cos (t) ^ 

v -sin(O y 


"sin(0" 

v cos(O y 


This is still a complex solution to a real differential equation. But were we 
to plug this expression into the equation Y = AY, the real part would remain 
real and the imaginary part would remain imaginary. That means that the 
real part of this complex solution, 


Y real (t)- 


( cos(f) ^ 
-sin(f) y 


and the imaginary part, 


Y imag (t) = 


"sin(0" 

v cos(O y 


are both solutions of this differential equation. And notice that the real part 
starts out at the point (1,0) whereas the imaginary part begins at (0, -1). So 
these are independent solutions. It is interesting that using just one of the 
eigenvalues gives us a pair of solutions that then generate all possible 
solutions (below) exactly as we saw earlier. 


Y(t) = k i 


r cos(£) ^ 

v -sin(0 y 


+ 


'sin(0' 

v c°s(0 y 


Clearly, when k { = 0, our solutions all lie on circles centered at the origin. In 
fact, that’s true for any values of k\ and k 2 (as long as they are both not 
equal to zero). So our phase plane here is seen in Figure 13.1. 
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Figure 13.1 


V 



This type of equilibrium is center. In general, if the eigenvalues are ±ia, the 
solutions may lie on circles or ellipses. 

When we compute the roots of a characteristic equation, we may also find 
complex roots of the form a + ib and a - ib where, unlike the previous case, 
a may be nonzero. In this case, the exact same procedure as above yields a 
complex eigenvector V for one of the eigenvalues. This then yields, by 
Euler’s formula, solutions of the form 

e (a+ib) ‘V = e a, (cos(bt) + ism(bt)) V. 

If we then break this vector solution into its real and imaginary parts, we 
find 2 real solutions to the linear system, and it turns out that these solutions 
are independent. So we find the general solution just as before. 

The difference here is that our new solutions always involve the exponential 
term e (at \ If a > 0, then this term tends to infinity as time goes on, whereas 
if a < 0, this term goes to zero as time goes on. As the remaining terms in 
the solution involve sin(Z^) and cos(Z?f), it follows that if a > 0, solutions 
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Lecture 13: Visualizing Complex and Zero Eigenvalues 


spiral away from the origin (the equilibrium is a spiral source), but if a < 0, 
solutions spiral in toward the origin (a spiral sink). 

For example, for the linear system 


we can easily check that the eigenvalues are -0.1 + i and -0.1 - i. Similar 
computations to the above then yield the general solution 


7(0 = v -0 ' 1 ' 


r cos(0 ' 
v -sin(0 y 


+ k 2 e 


/ sin(0' 
v cos(0 y 


The corresponding phase plane is the below. 



130 















If we consider instead the linear system 


7’ = 


0.1 1 
-1 0.1 


the eigenvalues are 0.1 + i and 0.1 - z, and the general solution becomes 


Y(t ) = V° 


cos(7) 

-sin(7) 


+ £ 2 e u 


^ sin(^) ^ 
cos(^) 


Solutions now spiral away from the origin, as seen below. 


Figure 13.3 


V 



So far we have figured out how to handle the major types of linear systems, 
sinks, saddles, sources, and centers. But there are some other special cases, 
such as when we have zero eigenvalues and real and repeated eigenvalues. 


If we have just one eigenvalue that is equal to zero (and so the other is 
either positive or negative), we know more or less how to solve the system. 
We find our eigenvector V corresponding to the eigenvalue 0, so we then 
have a solution of the form Y{t) = k\e 0t V. Our solution is just the constant 
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Lecture 13: Visualizing Complex and Zero Eigenvalues 


k\ V. So all points that lie on the straight line passing through the origin and 
the point V in the plane are constant or equilibrium solutions. Then the sign 
of the other eigenvalue a determines whether the remaining solutions tend 
toward or away from these equilibria. 

As an example, consider the system 


Check on your own that the eigenvalues are 0 and -1 and the general 
solution is 


Y(t) - k x e 


/A 


vV 


+ k o 




vO , 


So we have a line of equilibrium points along the x-axis, and all other 
solutions tend toward these equilibria. The phase plane is the below. 


Figure 13.4 


V 


^ W\t ^4 ^ Vt 

V ^ ^ ^ - 

^4^^4 ^4 iV ^ 

^4 ^ ^4 ^l\ 

* * 

'i ^ 'i 
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As we’ll see in the next lecture when we describe the trace-determinant 
plane, linear systems often undergo bifurcations when there is a family of 
such systems that has either a zero eigenvalue or a center. For example, 
consider the family of differential equations 


The eigenvalues here are -1 and A. So, whenever A < 0, we have a sink; 
whenever A > 0, we have a saddle. A bifurcation occurs when A = 0, and the 
system has a zero eigenvalue. 


Important Term 


spiral source: An equilibrium solution of a system of differential equations 
for which all nearby solutions spiral away from it. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 3.2-3.3. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 3.1. 
Kolman and Hill, Elementary Linear Algebra, chap. 7. 

Lay, Linear Algebra, chap. 5. 

Roberts, Ordinary Differential Equations, chap. 7. 

Strang, Linear Algebra and Its Applications, chap. 5. 

Strogatz, Nonlinear Dynamics and Chaos, chaps. 5.1-5.2. 
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Lecture 13: Visualizing Complex and Zero Eigenvalues 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Linear Phase Portraits. 


Problems 


1. a. For the system of differential equations 

( 2 5 ^ 

Y'= Y, 

1-1 -V 

first compute the characteristic equation. 

b. What are the eigenvalues of this matrix? 

c. Compute the asssociated eigenvectors. 


2. a. Compute the eigenvalues associated to 


7 ' = 


^0 


n 

-b 


Y. 


b. Find the associated eigenvectors for this system. 


3. Rewrite the second-order differential equation for the mass-spring 
system as a linear system and find the eigenvalues. Do you see any 
connection with our earlier method of solving this system? 
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4. a. 


What are the eigenvalues and eigenvectors of the matrix 


/ 


V 


a 

-b 


ZA 


? 


b. Determine the types of the equilibrium points for the linear system 
corresponding to the above matrix, and sketch the regions in the 
a-b plane where this matrix has different types of eigenvalues. 


5. a. Find the general solution of the system 


7 ' = 


' 1 3 > 

v V2 3V2 y 


Y 


and sketch the phase plane. 

b. Find the solution to this system satisfying the initial condition 

m = (i,o). 


Exploration 


Consider all possible 2x2 matrices of the form 


A = 


V c 


b\ 
d J 


Determine the set of ab-, o, and J-values for which this matrix has either 
a zero eigenvalue or pure imaginary eigenvalues (i.e., eigenvalues of the 
form ±p for some p ^ 0). Express these sets in terms of the trace and the 
determinant of A. 
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Summarizing All Possible Linear Solutions 

Lecture 14 


W e have solved almost every type of linear system. The only 
remaining case is that of repeated eigenvalues. The easy case in 
this situation is when our system assumes the form 


r f = 


a 0 
0 a 


Y . 


We see immediately that the eigenvalues are both given by a , and any 
nonzero vector is an eigenvector since 


'a (T 

x> 


f x l 

1° a ) 


= a 






Therefore all solutions are straight line solutions (assuming a is nonzero) 
either tending to the origin (if a < 0) or tending away (if a > 0). 

Here is the phase plane when a = 1. 


Figure 14.1 


v 
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This case, however, is unusual. Most often we find that the eigenvectors lie 
along a single straight line. To find the other solutions, consider the system 
given by 


It is easy to check that the eigenvalues here are -1 (repeated), and there is a 
single eigenvector given by (1, -1). So we have one straight line solution 
given by 


Y l (t) = k l e~‘ 


fl 

v-1 


\ 




To get the other solutions, recall that we have actually seen this example 
before; it is the linear system corresponding to the mass-spring system 
determined by the equation y" + 2y f +y = 0. We saw earlier that this was the 
critically damped harmonic oscillator whose general solution was given by 
y(t) = k\e f + k 2 te~ t . As a system, our solution is then given by 


Y(t) = 


r y(tf 


k { e ‘ + k 2 te ‘ ' 



K -k ie -‘ -k 2 te~ l +k 2 e~‘ 2 


which we can rewrite as 


Y(t) = k x e 


v-ly 


+ k 2 te 




v-l/ 


-\-k 2 £ 


vly 


Note that the eigenvector (1,-1) appears twice in this expression. But there 
is another term involving 
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Lecture 14: Summarizing All Possible Linear Solutions 


We easily check that the vector (0, 1) is actually the solution of the equation 


{A-XI) 


fX > 


rn 



-b 


where X is our known repeated eigenvalue. In fact, this is the prescription 
for finding the general solution for a linear system with repeated eigenvalue 
X and a single eigenvector. So the general solution is 

Y{t) = k x e M V + k 2 t e At V + k 2 e At W 


where we have 


(A-XI)W = V. 

The phase plane for the above system now has only one straight line 
solution, while all other solutions also tend to the origin since the terms e 1 
and te 1 both tend to zero as time goes on. 


Figure 14.2 


V 
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Indeed, as seen above, all other solutions tend to the origin in a direction 
tangential to the straight line solutions. The reason for this is that each of the 3 
terms in the solution involves either e~ f or te \ By calculus, as t tends to 
infinity, e 1 tends to zero much more quickly than te \ forcing all solutions to 
come into the origin tangent to the line containing the eigenvector. 

Now let’s summarize the entire situation for linear systems by painting the 
picture of the trace-determinant plane. For a linear system of the form 

fa b\ 

Y'= Y , 

\c d) 

we have seen that the characteristic equation is given by 
A 2 -TA + DA = 0. 


Here the trace T is given by a + d, and the determinant D is ad - be. Then 
the eigenvalues are given by 

r±Vr 2 -4D 

2 

Using this formula, we can summarize the situation for linear systems in 
writing as follows. 

1. The system has complex eigenvalues if T 2 - 4D < 0. So we have 

a. a spiral sink if T < 0, 

b. a spiral source if T > 0, and 

c. a center if T= 0. 
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2. The system has real, distinct eigenvalues if T 2 - 4D > 0. So we have 

a. a saddle if D < 0, 

b. a real sink if D > 0 and T > 0, and 

c. a real source if D > 0 and T < 0. 

3. The system has a single zero eigenvalue if D = 0 and T ^ 0. 

4. The system has repeated eigenvalues if T 2 - 4D = 0. 

We can also summarize this pictorially by drawing the trace-determinant 
plane. In this figure, the horizontal axis is the trace axis and the vertical axis is 
the determinant axis. Each matrix then corresponds to a point in this plane 
given by (7, D). In the figure below, the different regions then correspond to 
places where our linear system has different phase planes. The parabola given 
by T 2 - 4D = 0 represents the matrices that have repeated eigenvalues; the 
trace axis D = 0 is where we have a zero eigenvalue; and the positive 
determinant axis is where we have a center. The regions between these 3 
curves are where the other types of phase planes occur. 


Figure 14.3 
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This figure also represents the bifurcation diagram for linear systems. For if 
we have a family of linear systems that depends on a parameter k , then the 
corresponding curve that this family produces in the trace-determinant plane 
shows how and when bifurcations occur. For example, consider the family 
of linear systems given by 


The trace here is -1, and the determinant is k. So this family of systems lies 
along the vertical line given by trace = -1 and the determinant increasing. 
We see that this family starts out as a saddle, and then crosses the zero 
eigenvalue line and becomes a real sink. We eventually find repeated 
eigenvalues followed by a spiral sink. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps 3.5 and 3.7. 
Hirsch, Smale, and Devaney, Differential Equations , chaps. 3.3 and 4.1. 
Strang, Linear Algebra and Its Applications , chap. 5. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 5.2. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Linear Phase Portraits. 


141 








Lecture 14: Summarizing All Possible Linear Solutions 


Problems 


1. a. Sketch the direction field for the system 


Y' = 


r 2 

vO 


(T 

2, 


Y. 


b. What are the eigenvalues for this matrix? 

c. Find all possible eigenvectors for this matrix. 

d. What is the general solution of this system? 

e. How do typical solutions of this system behave? 


2. Find the general solution of the linear system 


1 (P 

i -b 


and sketch the phase plane. 


3. Find the general solution of the linear system 


Y' = 


( 0 


0" 

0, 


Y 


and sketch the phase plane. 
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4. a. 


Sketch the curve in the trace-determinant plane that the following 
family of linear systems moves along as a changes. 


7’ = 


OJ 


Y 


b. Find all a-values where this system undergoes a bifurcation (i.e., 
the type of the equilibrium point at the origin changes). 


5. Repeat problems 4a and 4b for the system 


Y' = 



2 ^ 


Y. 


Exploration 


Consider the 3-parameter family of linear systems given by 


7 ? = 




b > 
0 , 


7. 


First fix a > 0 and describe the analogue of the trace determinant plane (i.e., 
sketch the regions in the b-c planes where this system has different types of 
equilibria. Then repeat this question for a = 0 and a < 0. Put all of this 
information together to create a 3-dimensional model that describes all 
possible behaviors of this system. 
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Nonlinear Systems Viewed Globally—Nullclines 

Lecture 15 


W e now move on to nonlinear systems of ordinary differential 
equations. In this lecture we introduce one of the main 
techniques that we use to understand the behavior of solutions 
of nonlinear systems. This tool, called nullclines, allows us to partition the 
phase plane into certain regions where we know the approximate direction 
of the vector field. 

Given a nonlinear system x' = F(x , y) and y' = G(x, y ), the x-nullcline will 
be the set of points at which x' = 0—that is, F(x, y) = 0. At these points, 
the vector field points vertically, either up or down. Similarly, the 
y-nullcline is the set of points where / = 0, so G(x, y) vanishes. At these 
points, the vector field points horizontally, either to the left or the right. 
Note that the equilibrium points are then given by the points of intersection 
of the x- and y-nullclines. 

Most importantly, in the regions between the nullclines, assuming our 
vector field is continuous, the vector field must point in 1 of 4 directions: 
northeast, northwest, southeast, or southwest. This very often allows us to 
get a good handle on the qualitative behavior of nonlinear systems. 

Let’s start with a simple nonlinear system. 

x' = y — x 

y=i-y 

The x-nullcline is the parabola y = x 2 . The vector field is vertical along this 
parabola. It points to the right above the parabola since x' > 0 in this region, 
and it points to the left below the parabola. The y-nullcline is the horizontal 
line y = 1, so we immediately see equilibrium points at (-1, 1) and (1, 1). 
Moreover, the vector field is tangent to this nullcline. So we have solutions 
that move right along this line in the region above the parabola and to the 
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left on the portions below the parabola. We also have that y' < 0 above the 
lin Qy= 1 and thaty > 0 below this line. Therefore we know the direction of 
solutions in all the regions between the nullclines. 

If a solution starts in the region above the parabola but below the line y= 1, 
it must proceed up and to the right. It cannot cross the line y= 1, since we 
have a solution there and solutions cannot cross by the existence and 
uniqueness theorem. It also cannot hit the parabola because to do so, it 
must turn and become vertical. But the vector field always points upward 
here, so this is not possible. So solutions in this region must tend to the 
equilibrium point at (1, 1). The same thing happens in the region to the right 
of the parabola and above y =1; all these solutions tend to the same 
equilibrium point. In the left region below the parabola but above y = 1, 
solutions must tend to infinity. 

In the region above the parabola and y = 1, solutions have 4 choices: They 
can tend directly to (1, 1). They can cross the right branch of the parabola 
but then enter the region where all solutions tend to (1, 0). They can hit the 
left branch of the parabola, in which case they tend to infinity. Or they can 
tend to the equilibrium point at (-1, 1). In the lower region, solutions 
similarly have 4 possible choices, but we know more or less what happens 
to all solutions. It appears that (1, 1) is a sink and (-1, 1) is a saddle. But the 
latter is not clear: Can we have more than a single curve tending to that 
equilibrium point? We’ll see how to handle this in the next lecture. 

Let’s go back to the predator-prey system 

x' = x(l ~ x) - xy 
y' = ~y + axy. 

Here we first choose 0 < a < 1, so the rate of predation is small. Recall that 
the computer showed us that all solutions tended to the equilibrium point at 

x = l,y = 0. 
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Figure 15.1 
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The x-nullcline will be the set of points in the phase plane where the 
x-component of the vector field vanishes. So the x-nullcline is the set of points 
where the vector field points vertically. In this example, the x-nullcline is 
given by x(l - x) - xy = 0, so that xy = x( 1 - x). Thus the x-nullcline consists 
of the pair of straight lines x = 0 and y = 1 — x. Similarly, the y-nullcline is the 
set of points where / = 0 (i.e., where the vector field points horizontally). In 
this case, the y-nullclines are the lines y = 0 and x = Ha. 


Figure 15.2 

y 
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In this figure, the solid lines are the x-nullclines and the y-nullclines are 
dashed. As before, the places where an x- and a y-nullcline intersect are the 
equilibrium points for the system. So in this case, the equilibria are (0, 0) 
and (1, 0). Note also that the y-nullcline x = 1 /a and the x-nullcline y = 1 —x 
meet at a point where the y-coordinate is negative. This means that this is 
not an equilibrium point for the predator-prey system, since both x and y are 
nonnegative. 

Most important, at all points in one of the regions between the nullclines, 
the vector field always points in 1 of 4 directions: northeast, northwest, 
southeast, or southwest. And we can determine this direction by simply 
computing the vector field at a single point in the region. For example, at 
(1, 1), our vector field is given by x' = —1, y' = — 1 + a < 0; so the vector 
field points down and to the left at all points between the lines y= 1 - x and 
x = Ha. In the triangle bounded by y = 1 - x and the x- and y-axes, the 
vector field points down and to the right. And in the region to the right of 
x = 1/a, the vector field points up and to the left. 


Figure 15.3 
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As a consequence, we completely understand what happens to all solutions 
of this system. Solutions that start to the right of x = Ha must move up 
and to the right until they hit the line x = Ha. Then they move down and 
to the left. So there are only 2 possibilities for what happens to these 
solutions: Either the solution tends to the equilibrium point at (1, 0) or it 
crosses the line y = 1 — x. (They can never cross the x- or y-axes by the 
uniqueness theorem.) 

In the final case, in the triangular region, solutions move down and to the 
right. Since these solutions can never cross the x-axis, they must also tend to 
the equilibrium point at (1, 0). We therefore know that all solutions of this 
predator-prey system necessarily tend to the equilibrium point at (1, 0), 
assuming that our initial populations are both positive. Therefore we know 
that when a < 1, the predator population goes extinct. 

Nullclines do not always give the complete picture of the phase plane. For 
example, if a > 1 in our predator-prey system, we get the following 
nullcline picture. Here the nullclines x = Ha and y = 1 - x do cross at a new 
equilibrium point ( Ha , (a - 1 )/a). 


Figure 15.4 
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We do know that every solution starting at a point with x, y > 0 must 
eventually enter the region x < 1 la. But does this solution tend to the new 
equilibrium point? If so, does the solution spiral into this equilibrium point, 
or does it not spiral? The computer seems to indicate that both types of 
behavior can happen. These figures below show the phase plane when 
a = 1.2 (left) and when a = 2 (right). 


Figure 15.5 



v 



In addition, could the solution move around the equilibrium point in a 
periodic motion and never come to rest? This kind of local behavior near 
equilibria will be the subject of the next lecture. 

As another example of nullcline analysis, consider the nonlinear system 

/=*-/• 

The x-nullcline for this system is the parabola y = x 2 , which opens upward; 
the y-nullcline is the parabola x = y 2 , which opens to the right. These 
nullclines meet at the 2 equilibria: (0, 0) and (1, 1). In the 5 regions 
separated by the nullclines, we compute that the vector field points 
as follows. 
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Figure 15.6 



So (1, 1) appears to be a sink and (0, 0) a saddle. We don’t know this for 
sure, but the next lecture will give us the local tool to make this assessment. 
The actual phase plane for this system is below. 


Figure 15.7 
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Important Term 


existence and uniqueness theorem: This theorem says that, if the right- 
hand side of the differential equation is nice (basically, it is continuously 
differentiable in all variables), then we know we have a solution to that 
equation that passes through any given initial value, and, more importantly, 
that solution is unique. We cannot have two solutions that pass through the 
same point. This theorem also holds for systems of differential equations 
whenever the right side for all of those equations is nice. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 5.2. 
Edelstein-Keshet, Mathematical Models in Biology , chap. 6.2. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 9.1. 
Strogatz, Nonlinear Dynamics and Chaos , chap. 6.1. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Linear Phase Portraits. 


Problems 


1. a. Find the x- and y-nullclines for the system 

x' =y 
/ = 2x. 
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b. Sketch the regions where the vector field points in 
different directions. 

c. What are the equilibrium points for this system? 

d. Sketch the direction field and phase plane for this system. 

e. How many solutions tend to the equilibrium points? 


2. Sketch the nullclines and phase plane in the region x, y > 0 for the 
following system. 

x' = x(—x ~ 3y + 150) 
y'=y(-2x~y+ 100) 

3. Sketch the nullclines and phase plane in the region x, y > 0 for the 
following system. 

x' = x(100 - x - 2 y) 
y' =y(150 - x - 6y ) 

4. Using nullclines, determine the phase plane for the following system. 

x' = x(l - x) 

/=*-/ 
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5. a. Consider the predator-prey system 
x' = x(l — x) ~ Axy 
y'=y( i ~y)+xy, 

where A is a positive parameter. Using nullclines, determine the 
values of A where different outcomes occur. In particular, 
determine at which parameters this system undergoes a bifurcation. 

b. What changes will there be in the predator-prey system if we 
introduce a new parameter of B instead of A? 

x' = x(l — x) — xy 
y'=y( 1 ~y) + Bxy 


Exploration 


Use the nullclines to investigate the phase plane for the following system 
that depends on a parameter A. 

x' = x 2 ~\ 
y = ~xy + A(x 2 - 1) 

At what values of A do you see a bifurcation? Explain what happens to 
solutions before and after the bifurcation. 
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Nonlinear Systems near Equilibria—Linearization 

Lecture 16 


Recall that for the predator-prey system in the last lecture 

x' — x(l - x) - xy 
y' = -y + axy 

we had an equilibrium point at x = 1 /a,y = (a - 1 )la when a > 1. It appeared 
from the computer illustrations that sometimes solutions spiraled in to this 
point, but for other a-values this did not occur. To understand better what 
happens near equilibria in nonlinear systems, we use a procedure called 
linearization. This often allows us to use a linear system to figure out what 
is happening locally near a nonlinear equilibrium point. 

To explain the process of linearization, let’s first consider a first-order 
equation y' =fly), and let’s assume we have an equilibrium point aty = 0. 
Recall from calculus that we can expand the function fly) as a Taylor series 
about y = 0. That is, we can write it as follows. 

Ry) = /(0)+ fXVy+^-y 2 +- 

The first term in this series is 0, since we have an equilibrium point aty = 0. 
Also, if y is very close to 0, they 2 term (the third term in the series) is much 
smaller than the term involving just y (the second term). Similarly, the 
terms involving y n for n > 2 are even smaller. This tells us that we can 
approximate the right-hand side ofy' =fly) by just using the second term in 
the Taylor series, namely/'(0)y. Therefore, very close toy = 0, solutions of 
our nonlinear ODE y' = fly) should be close to solutions of the ODE 
y = f'($$)y, which is just a linear ODE that we know how to solve. 

If our equilibrium point is located at a nonzero value, say y 0 , then a simple 
change of variables says that the same process is true: Solutions of the 
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nonlinear ODE resemble the solutions of the linear equation;;' =f'(yo)y , at 
least in a very small neighborhood of y = y 0 . There is one slight problem 
here. If / '(y 0 ) happens to equal 0, then the linear system reduces to just 
y = 0, for which all solutions are equilibria. This certainly will not be the 
case for the corresponding nonlinear system. 

For a nonlinear system 

x’ = F(x,y) 
y' = G(x,y) 

with an equilibrium point at (x 0 , y 0 ), we can perform a similar linearization, 
only now the corresponding linear approximation involves the partial 
derivatives of both F and G. That is, our nonlinear system has solutions near 
(x 0 , y Q ) that resemble the solutions of the following linear system. 

dF dF 

*' =— (x 0 ,y 0 )x+—(x 0 ,y 0 )y 
ox dy 


, dG . dG 

y =^~{x 0 ,y 0 )x +—{x 0 ,y 0 )y 
dx oy 


Note that this is a linear system of ODEs. Again, what we have done here is 
just drop all the higher-order terms in the Taylor series expansions of both F 
and G. As above, solutions of our nonlinear system near the equilibrium 
point (x 0 , y 0 ) should resemble the solutions of the linear system 

Y' = JY, 


where J is the matrix 


f dF dF y 

-T-OWo) -^-(Wo) 


dx 

dG 


dy 


a (^>^0) ir-(wo) 

ox oy 
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The matrix J is called the Jacobian matrix for the nonlinear system. 


As an example, consider the system below. 

x' = x — xy 2 
y' = ~y + xy 

There is no way to solve this system explicitly, but we certainly have an 
equilibrium point at the origin. The Jacobian matrix at an arbitrary point is 

(+2 \ 

J= 1 ~y -2 xy 

v y -i+x. 

At the origin we get 

(1 0 \ 

J = 

to -lj 

so our nonlinear system is close to T = JY. Since J has eigenvalues 1 and 
-1, our nonlinear system has a saddle at the origin. A glance at the phase 
plane shows that there is a lot more going on in this system. 


Figure 16.1 
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However, if we zoom in at the origin, we do see a part of the phase plane 
that looks like a saddle. 

Figure 16.2 


v 



x 


05 


As in the first-order case, there are times when linearization does not work. 
In particular, if the Jacobian matrix associated to our equilibrium point has 
imaginary or zero eigenvalues, then the addition of the higher-order terms 
will usually greatly alter the linear phase plane. For example, the system 

x' =y — x(x 2 + y 2 ) 
y’ = -x-y(x 2 + y 2 ) 

has an equilibrium point at the origin. The corresponding Jacobian matrix is 
given by 



J = 


and the eigenvalues are easily seen to be the imaginary numbers i and ~i. So 
the linearized system has a center at the origin. But this is not true of the 
nonlinear system (Figure 16.3), where we see a spiral sink. 
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Figure 16.3 


v 



If we look instead at the system 

x' = y + x(x 2 + y 2 ) 
y' = -x+y(x 2 + y 2 ) 

then again the linearized system has a center at the origin. But the same 
arguments as above show that this equilibrium point is now a source. 

Let’s return to the predator-prey system 

x' = x(l —x)—xy = F(x, y ) 
y = —y + axy — G(x, y). 


As we saw earlier, we have an equilibrium point at (1 la, (a - 1 )/a) when 
a > 1. The Jacobian matrix at an arbitrary point is 


J = 


^1-2 x + y 
v a y 


-x 

-1 + ax 


\ 

) 
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At our equilibrium point, this matrix becomes 


f-l/a -1 / a\ 


so the solutions near our equilibrium point should resemble those of the 
linearized system Y' = JY. The characteristic equation for the linearized 
system is 

T + /L / cl + (cl — 1) / a = 0 , 
so the eigenvalues are 

X 0 (-l±Vl + 4a(l-«)). 

There are then 2 different cases. The first occurs when the term inside the 
square root is negative. In this case, the eigenvalues are complex with 
negative real parts, so the equilibrium is a spiral sink. If the term 1 + 4a(l - 
a) > 0, we certainly know that this term is less than 1. This is because a > 1, 
which makes Aa{\ - a) negative. Our eigenvalues are both real and 
negative, so in this case we have a real sink. An easy computation shows 

that the first case occurs when a > 2 + 2^2 , and the second case occurs if 
\< a <2 + 2\[2 . This is what we saw in the previous lecture. 


Important Terms 


Jacobian matrix: A matrix made up of the various partial derivatives of the 
right-hand side of a system evaluated at an equilibrium point. The 
corresponding linear system given by this matrix often has solutions that 
resemble the solutions of the nonlinear system, at least close to the 
equilibrium point. 
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Taylor series: Method of expanding a function into an infinite series 
of increasingly higher-order derivatives of that function. This is used, 
for example, when we approximate a differential equation via the technique 
of linearization. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 5.1. 
Edelstein-Keshet, Mathematical Models in Biology , chaps. 4.7-4.8. 
Hirsch, Smale, and Devaney, Differential Equations , chaps. 8.1-8.3. 
Strogatz, Nonlinear Dynamics and Chaos , chap. 6.3. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , HPG Linearizer, HPG System 
Solver, Vander Pol. 


Problems 


1. Compute the partial derivatives with respect to x and y of the function 

F(x,y) = (x 2 y+x i + y 3 ). 

2. What is the Jacobian matrix for a linear system 

x' = ax + by 
y = ex + dyl 









3. a. Find all equilibrium points for the system 

x' = x 2 + y 
y=x. 

b. Compute the Jacobian matrix at each equilibrium point. 

c. What is the type of these equilibrium points? 


4. Find and determine the type of the equilibrium points for the Vander 
Pol equation given by 

x’ =y 

y' = -x+( 1 ~x 2 )y, 

and then use the computer to see what else is going on for this system. 


5. Give an example of a nonlinear system that has an equilibrium point for 
which the linearized system has a zero eigenvalue, but the nonlinear 
system has only one equilibrium point. 


6. For the system 

x'=y~l 

y’=y-x 2 , 

determine the types of each equilibrium point. 
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7. 


For the system 


x' = sin(x) 
y = cos(y), 

find all equilibrium points and use linearization together with the 
nullclines to sketch the phase plane. 


8. For the system 

x' =y~x 2 +A 

y=y, 

determine the types of each equilibrium point and also the ^-values 
where a bifurcation occurs. 


Exploration 


Consider the system of differential equations 

x' = x/2 - y - (x 3 + y 2 x)/2 
y = x + y/2 ~ (y 3 + xy)/ 2. 

What is the type of the equilibrium point at the origin? Using the computer, 
find out what happens far away from the origin. To understand this 
behavior, change coordinates to polar coordinates. Then you will be able to 
solve this system explicitly. 
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Bifurcations in a Competing Species Model 

Lecture 17 


L et’s combine the previous topics to investigate what happens when 
we have 2 species that compete with each other for the same 
resources. Let x and y denote the population of these 2 species. We 
first assume that if the other species is absent, the population of the given 
species obeys the limited-growth population model. So when y = 0, we will 
assume that the x-population is given by x' = x(l - x/M). When x = 0, we 
havey' = y( 1 ~y/N). So our phase plane thus far is the below. 


Figure 17.1 



Our second assumption is that if the y population increases, it should have a 
negative effect on the x-population, and vice versa. So if y increases, x' 
should decrease. One way to model this is by adding the term —(a/M)xy to 
the differential equation for x and ~(b/N)xy to the y-equation. So our 
competing species system is the below. 

dx 

— = x(l - x / M) -{a! M)xy 

— = y(l - y / N) — (b/ N)xy 

dt 


We have 4 parameters for this system: a , b , M, and N. The quantities a and 
b measure how competitive the 2 species are. For simplicity, let us assume 
that M = N = 400; the equations are below. (Note: There is nothing special 
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about the number 400; it is simply the default value of M and N in the 
software we use.) 

dx 

— = x(l-;t/400-ay/400) 
dt 

—= y(l — y 1 400 — bxl 400) 
dt 

Now let’s view the phase plane and nullclines. First assume that a, b > 1 
(so our species are very competitive). The x-nullclines lie on the y-axis 
(where we know what happens) and along the straight line 1 - (x/400) - 
ay 1400 = 0, or y = 400/a - x/a. Similarly, the y-nullclines are given by the 
x-axis and the line y = 400 - bx. Note that the 2 nullclines that do not lie on 
the axes meet at the point x = 400(1 - a)/( 1 - ab ), y = 400(1- b( 1 - a)/ 
(1 - ab)). This is another equilibrium point for our system since a, b > 0. 
We call this the coexistence equilibrium point. 

At the point (400, 400), we compute that x' = -400a and / = -400 b, so the 
vector field points down and to the left. Similarly, at (1, 1) we have x' > 0 
and y > 0, so the vector field points up and to the right. We similarly 
compute the value of the vector field in the other regions bounded by the 
nullclines to obtain the following figure. 


Figure 17.2 
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It appears that the coexistence equilibrium point is a saddle. Of course, we 
do not know this for sure from the picture. However, linearization at this 
point should give us the answer. For simplicity, assume a = b = 2; so our 
equations are 

x' = x- x 2 /400 - xyl .200 
y' = y~ y 2 /400 - xy/ 200 . 

The coexistence equilibrium point is found by solving the 2 equations 

1 - x/400 - y/200 = 0 
1 - y/400 - x/200 = 0 

simultaneously. A little algebra yields x = y = 400/3. Then the Jacobian 
matrix is 

r \ — x! 200-y / 200 -x/200 > 

v -y/ 200 1-x / 200-y / 200J ’ 

and at (400/3, 400/3) this matrix is 

'-1/3 -2/3 > 
v -2/3 -1/3/ 

The determinant of this matrix is -1/3. Since the determinant is negative, 
we know that the point in the trace-determinant plane that corresponds to 
this matrix lies below the trace axis. So we are indeed in the region where 
we have a saddle. The phase plane for a, b> 1 is shown in Figure 17.3. 
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Figure 17.3 
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Any solution that starts at a point (x 0 , yo) with x 0 and y 0 positive (except for 
the special solutions that tend to the saddle) tends to either (400, 0) or 
(0, 400). That is, since both species are very competitive, one will almost 
surely prevail and the other will go extinct. Which species prevails depends 
on the initial populations of the 2 species. 

Now suppose that species y is much less competitive, so the constant a 
decreases. Recall that one of the x-nullclines is given by y = 400 la - x/a. So 
the y-intercept of this line is 400 la. Note that this y-intercept agrees with the 
y-intercept of the y-nullcline given by y = 400 - bx, which is 400 when a is 
equal to 1. So we see a bifurcation at a = 1 since the x- and y-nullclines that 
do not lie on the axes now meet at the point (0, 400). That is, the 
coexistence equilibrium point merges with the equilibrium point at (0, 400). 
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Figure 17.4 
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One can check that the Jacobian matrix at the equilibrium point (0, 400) is 
given by 

(0 0 ^| 


This matrix has eigenvalues -1 and 0, so linearization does not give us an 
accurate picture of the phase plane. However, the configuration of the 
nullclines shows that all solutions that do not lie on the y-axis now tend to 
the equilibrium point on the x-axis at (400, 0). That is, the species y now 
goes extinct no matter what the initial point (x 0 , yo) is (as long as x 0 
is positive). 

When a < 1, the nullclines again indicate that all solutions tend to (400, 
0)—the weaker species y remains extinct. Similarly, if b becomes less than 
1 while a stays larger than 1, we see that the x-population goes extinct. 
Below is the nhase nlane for a < 1. b > 1. 
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Lecture 17: Bifurcations in a Competing Species Model 


Figure 17.5 
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Next, let’s look at what happens if both a and b are small (< a , b < 1). Both of 
our 2 species are not that competitive. We see another bifurcation as b moves 
below 1 while a stays at that level. A new coexistence equilibrium point is 
bom. The nullclines, however, indicate that this equilibrium is now a real 
sink. Indeed, the lack of strong competition leads to coexistence no matter 
where the 2 populations start out. The below is the phase plane when a, b< 1. 
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Finally, we can record all of this information in the analogue of the trace- 
determinant plane or the bifurcation diagram by plotting the different 
regions in the a-b plane where we have different behaviors. This also allows 
us to see where bifurcations occur. 


Figure 17.7 



Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 5.1-5.2. 
Edelstein-Keshet, Mathematical Models in Biology , chap. 6.3. 

Hirsch, Smale, and Devaney, Differential Equations , chap. 11.3. 
Strogatz, Nonlinear Dynamics and Chaos , chap. 6.4. 
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Lecture 17: Bifurcations in a Competing Species Model 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Competing Species, HPG 
System Solver. 


Problems 


1. a. Consider a simpler version of the competing species model where 

we do not include overcrowding: 

x' = x - xy 
y'=y-xy. 

Here both x and y are non-negative. First, find all equilibrium 
points of this system. 

b. Compute the Jacobian matrix at each equilibrium point. 

c. Determine the types of these equilibria. 

d. Determine the nullclines and sketch the regions where the vector 
field points in different directions. 

e. Sketch the phase plane for this system with x, y > 0. 

2. a. Verify that the coexistence equilibrium point in the competing 

species model is indeed a real sink when a, b= 1/2. 

b. Use linearization to check the types of the equilibrium points that 
lie on the x- andy-axes for all a- and 6-values. 







3. a. A model due to Beddington and May for the competitive behavior 
of whales and krill (a small type of shrimp) is 

x f = x(l —x)—xy 
y'=y( i ~y/x). 

Here y is the population of whales. The carrying capacity for the 
whale population is not constant; rather, it is given by x and 
therefore depends on the krill population. Find the equilibrium 
points for this system. 

b. Determine the types of the equilibria for the system above. 

c. Use nullclines and the computer to paint the picture of the phase 
plane for this system. 


Exploration 


Here is a modified competing species model where we now allow 
harvesting or immigration (the parameter h). 

x' = x( 1 — ax—y) 
y = y(b -x—y) + h 

First assume that h = 0 (no harvesting). Give a complete synopsis of the 
behavior of this system by plotting the regions in the a-b plane where 
different behaviors occur. Then repeat this for the parameters h > 0 and h < 
0. Try to envision the complete picture by viewing the regions in the a-b-h 
space where different outcomes occur. 
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Lecture 18: Limit Cycles and Oscillations in Chemistry 


Limit Cycles and Oscillations in Chemistry 

Lecture 18 

F or most of the systems of differential equations we have seen thus 
far, the most important solutions have usually been the equilibrium 
solutions. But there is another type of solution that often plays an 
important role in applications: limit cycles. These are solutions Y(t) that 
satisfy Y(t + A) = Y(t) for all values of t and have the additional property 
that they are isolated. By isolated, we mean that there are no nearby limit 
cycles; all nearby solutions either spiral in toward the limit cycle or spiral 
away from it. 

As an example of a limit cycle, consider the system 

x' = -y + x( 1 - (x 2 + y 2 )) 
y'=x+y(l-(x 2 +y)). 

Note that this vector field can be broken into a sum of 2 pieces: the 
linear system x f = ~y, y f = x; and the nonlinear system given by x' = 
x(l - (x 2 + y 2 )), y = y{ 1 - (x 2 + y 2 )). The linear system is a center with 
solutions traveling around circles centered at the origin. The nonlinear 
system is a vector field that points directly away from the origin inside the 
unit circle, points directly toward the origin outside the unit circle, and 
vanishes on the unit circle. So when we add these 2 vector fields, we find 
that there is one solution that travels periodically around the unit circle, 
while all other solutions (except for the equilibrium point at the origin) 
spiral toward this periodic solution. 

So the solution lying on the unit circle is a stable limit cycle, since nearby 
solutions move toward this solution. If we change the +x and +y to —x and 
—y in the second vector field, the solutions now spiral away, so we have an 
unstable limit cycle. 
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Figure 18.1 
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One of the principal tools for finding limit cycles is the Poincare-Bendixson 
theorem, which says the following: Suppose you have a ring-shaped region 
A in the plane where the vector field points into the interior of A along each 
of the 2 boundaries. Then provided there are no equilibrium points in A (and 
the vector field is sufficiently differentiable), there must be at least one limit 
cycle in A. In the figure below, we see a ring-shaped region A with limit 
cycle y contained inside. 


Figure 18.2 
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Lecture 18: Limit Cycles and Oscillations in Chemistry 


One of the most interesting discoveries in chemistry occurred in the 1950s, 
when the Russian biochemist Boris Belousov discovered that certain 
chemical reactions could oscillate back and forth between different states 
for long periods. Back then, most chemists thought that all chemical 
reactions settled directly to their equilibrium states. (For the interesting 
history, see the Winfree article in the Suggested Reading.) 

Since that time, many chemical reactions have been found to oscillate. 
We’ll look at a simpler model involving chlorine dioxide and iodide. The 
nonlinear system is 




where x and y are the concentrations of iodide and chlorine 
dioxide, respectively. There are 2 parameters here, A and B. We will 
simplify things by setting A = 10. A little algebra shows that there is a 
single equilibrium point at (2, 5). Linearizing at this equilibrium point gives 
the Jacobian matrix 

^ 7/5 —8/5 ^ 

v 85/5 -25/5/ 

The trace of the Jacobian matrix is T = -y-J B, and the determinant is 

D - 2 B. Note that T = 0 and D = 1 when B - 3.5. So in the 
trace-determinant plane, we see that the linearized equation has a center 
when B = 3.5. When B is a little larger than 3.5, we have T < 0, so the 
equilibrium point is a spiral sink; when B is slightly smaller than 3.5, the 
equilibrium point is a spiral source. So we have some sort of bifurcation 
going on when B = 3.5. What happens both mathematically and chemically? 
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For the whole picture, we plot the nullclines. The x-nullclines are given by 
x = 0 andy = 1 + x 2 , and the y-nullcline is given by 

_ (10-x)(l + x 2 ) 

7 " 4x 

Note that the y-nullcline has a vertical asymptote along the y-axis and 
crosses the x-axis only at x = 10. As usual, we compute the vector field in 
the regions between the nullclines to find what the phase plane looks like. 

It is clear in this graph that solutions spiral around the equilibrium point. 
But how do they do this? When B > 3.5, we have a spiral sink at (2, 5), and 
the computer shows that all solutions spiral into the equilibrium point. 


Figure 18.3 
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Lecture 18: Limit Cycles and Oscillations in Chemistry 


When the parameter B goes below 3.5, the equilibrium point changes to a 
spiral source—but something else happens. If we draw the rectangle 
bounded by the axes and the lines x = 10 and y = 101, we see that the vector 
field points into the rectangle R along all of the boundary points. Meanwhile 
we have the spiral source at (2, 5), so we can find a small disk D around this 
point so that the vector field points outside of D at all points on the 
boundary of D. So we have a ring-shaped region R- D on whose boundary 
the vector field points into the region R - D. By the Poincare-Bendixson 
theorem, there must be a limit cycle inside R — D. 

In fact, what has happened is that this system has under gone a Hopf 
bifurcation at B = 3.5. Mathematically, when ^ > 3.5, all solutions tend to 
the spiral sink. But as B passes through 3.5, not only does the spiral sink 
become a spiral source, but we also see the birth of a limit cycle. In fact, all 
solutions (except the equilibrium point) now tend to this limit cycle, so this 
is a stable limit cycle. 


Figure 18.4 
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Chemically, when B > 3.5, the reaction oscillates down to its equilibrium 
state. But when B> 3.5, the reaction oscillates back and forth forever. 


Important Terms 


Hopf bifurcation: A kind of bifurcation for which an equilibrium changes 
from a sink to a source (or vice versa) and, meanwhile, a periodic solution 
is bom. 

limit cycle: A periodic solution of a nonlinear system of differential 
equations for which no nearby solutions are also periodic. Compared with 
linear systems, where a system that has one periodic solution will always 
have infinitely many additional periodic solutions, the limit cycle of a 
nonlinear system is isolated. 


Suggested Reading 


Edelstein-Keshet, Mathematical Models in Biology , chap. 8.8. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 10.7. 
Strogatz, Nonlinear Dynamics and Chaos , chap. 8.3. 

Winfree, “The Prehistory of the Belousov-Zhabotinsky Reaction.” 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Chemical Oscillator, HPG 
System Solver. 
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Lecture 18: Limit Cycles and Oscillations in Chemistry 


Problems 


1. a. Does the system 

x'=x 2 +y 2 
y=x 2 + y 2 

have any limit cycles? Explain why or why not. 
b. What is the behavior of solutions of this system? 


2. Show that the system 

x f = x —y - x(x 2 + y 2 ) 
y’=x+y-y{x 2 +y 2 ) 

has a periodic solution on the unit circle. 


3. Does the system 

x' = x - x(x 2 + y 2 ) 
y' = x~y(x 2 +y 2 ) 

have a periodic solution on the unit circle? 


4. Check that the equilibrium point for the oscillating chemical reaction in 
the lecture summary is indeed (2, 5). 
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5. Describe the behavior of solutions of each of the following systems 
given in polar coordinates. 

a. r' = r-f,Q'= 1 

b. r' = ? - 1? + 2 r, 0' = 1 

c. r' = sin(r), 0' = -1 


6. Describe the bifurcations that occur in the system given by r' = ar - 
r 2 + r 3 , 0' = 1, where we assume that a > 0. 


Exploration 


A chemical reaction due to Schnakenberg is determined by the system 

x' =x 2 y ~x - 1/10 
y = ~x 2 y + a , 

where a is a parameter. Find the equilibrium points for this system, and 
compute the linearization. Do you observe any bifurcations? Can you show 
the existence of limit cycles for certain parameters? 
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All Sorts of Nonlinear Pendulums 

Lecture 19 


I n this lecture, we introduce several new types of methods used to 
understand nonlinear systems of differential equations. These include 
Hamiltonian and Lyapunov functions. In each case, we will concentrate 
on various nonlinear pendulum equations. 

The pendulum equations arise as follows. Suppose we have a light rod of length 
L with a ball (the bob) of mass m at one end. The other end of the rod is 
attached to a wall and is free to swing around in a circular motion. Our goal is to 
describe the motion of the ball. Suppose the position of the center of mass of the 
ball at time t is given by 0(f). Let’s say that 0 = 0 is the downward-pointing 
direction and that 0 increases in the counterclockwise direction. 

We will assume that there are 2 forces acting on the pendulum. The first is 
the constant force of gravity given by mg (where g « 9.8 m/s 2 ). Only the 
force tangent to the circle of motion affects the motion, so this force is 
-mgsin(0). We determine this with trigonometry. 


Figure 19.1 
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The position of the bob at time t is (Zsin(0(7)), -Zcos(0(7))). The speed of 
the bob is then the length of the velocity vector, which is given by 
(Zcos(0(7))0'(7)> Tsin(0(7))0'(O)- So the length of this vector is Z0'. 
Similarly, the component of the acceleration vector in the direction of the 
motion is L0". Our second force is then the force due to friction, which we 
assume as in the case of the mass-spring system is proportional to velocity. 
So this force is given by ~bLQ'(t ), where b is again called the 
damping constant. 

Newton’s second law, F = ma, then gives the second-order equation for 
the pendulum: 

9"+ — 9'+—sin(9) = 0. 

m L 


As a system, we get 


6' = v 


v f — 



fsin (9). 


Note that because of the sine term in the second equation, this is a nonlinear 
system of differential equations. 

Let’s first consider the case where there is no friction. This is called the 
ideal pendulum. We assume as usual that L = m = 1, so the system of 
equations becomes 

6' = v 

v’ = -gsin(0). 
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As always, we first find the equilibrium points. Clearly, they are given by v 
= 0 and sin(0) = 0. So our equilibria lie at v = 0 and 6 = nn (where n is any 
integer). When n is even, this equilibrium point corresponds to the 
pendulum hanging at rest in the downward position. When n is odd, the 
pendulum is at rest in its perfectly balanced upward position. There is a 
conserved quantity for this system of equations. Consider what physicists 
call the energy function: 

H(0,v) = \v 2 -gcos(d). 


In our solutions, both 0 and v are functions of t, so we can consider the 
energy function also to be a function of t. And then we can compute its 
derivative with respect to t. We find that 


dH 

dt 


vv g sin($) O' -vg sin(<9) - vg sin(<9) = 0 . 


Therefore H must be constant along any of the solution curves of the 
system. So we just plot the level curves of the function H (i.e., the curves 
given by H= constant). Then the solutions must lie along these level curves. 


Figure 19.2 
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We know that our vector field is everywhere tangent to these level curves, 
so our solutions must therefore look like the below. 


Figure 19.3 


V 



Besides our equilibrium solutions where the pendulum is either perfectly 
balanced in the upright position or hanging motionless in the downward 
position, we see 3 kinds of other solution curves labeled^, B , and C. 


Figure 19.4 



Close to the origin (or to any equilibrium point of the form (Inn, 0)), we 
see circular motions labeled A. These solutions correspond to the pendulum 
swinging back and forth periodically without ever crossing the upward 
position 0 = n. The second type of solution (labeled B) is where the 
pendulum moves continuously in the counterclockwise direction (or in the 
clockwise direction if v < 0). 
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Lecture 19: All Sorts of Nonlinear Pendulums 


The solutions labeled C tend to the upright equilibrium points as t tends to 
both ±oo. These upright equilibria are given by (±n, 0) . Using 

linearization, we can easily check that all of these points are saddle points. 
So this solution is what is called a separatrix (or saddle connection): It 
tends from one saddle point to another. These are the solutions that you 
would never see in practice. As time goes on, the pendulum tends to the 
upright equilibrium without ever passing through 0 = n. 

The ideal pendulum is an example of a Hamiltonian system, named for the 
Irish mathematician William Rowan Hamilton (1805-1865). A Hamiltonian 
system in the plane is a system of differential equations determined by a 
function H(x , y), which is called the Hamiltonian. The system is then given by 


x' = ^-(x,y) 
dy 


y' 


dH_ 

dx 


(x,y). 


H is a conserved quantity since, as before, if we differentiate H with respect 
to t, we find that 


dH dH dH dH dH dH ( dH\ ^ 

— =—jc ! (0+— y\t) = -+—-=o 

dt dx dy dx dy dy \ dx) 


We have seen another Hamiltonian system earlier, namely, the undamped 
mass-spring system: 

y = v 

V = ~ay. 
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Here the Hamiltonian is given by 


H{y,v)=\v 2 +±y 2 

since 


y 


_ dH 

~ V ~ dv 


and 


v' = -ay - - 


dH 

~sy' 


All level curves of H here are just ellipses or circles surrounding the origin 
(i.e., we have a center equilibrium point at the origin). 


Now let’s return to the original pendulum equation, but this time we will 
allow friction. Our second-order equation was 

6 "+ bO ’+ g sin($) = 0 , 

assuming L = m = 1 or, as a system, 

6' = v 

v’ = -6v-gsin(0). 
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This system is no longer Hamiltonian. But look at our former energy 
function from the ideal case 


H(6,v) = \v 2 -g cos(0). 


We compute 

- v-v'-gsin^)#' = v(-Z?v-gsin(6 > ) ) + gsin(<9) v = - bv 2 < 0 . 


This says that H(0(t),v(t)) is now a nonincreasing function along solution 

curves. So our solutions must now descend through the level curves of the 
function H. So we again know what happens to our solutions. Recall that 
the level curves of H looked like the following. 

Figure 19.5 


V 
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Now our solutions must look something like the below. 


Figure 19.6 


V 


0 



We see that now all solutions tend to the equilibrium point at which the 
pendulum hangs straight downward. A function like H is called a 
Lyapunov function, which is named for the Russian mathematician 
Aleksandr Lyapunov (1857-1918). 


You may ask, when is a system of the form 

X = F(x, y) 

Y = G(x, y) 

Hamiltonian? We would need F(x, y) = SH/dy and G(x, y) = -dH/dx. If such 
an H exists and is twice continuously differentiable, we would need to see 
the following. 

-dG _ d 2 H _ d 2 H _ dF 
dy dydx dxdy dx 
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Important Terms 


Lyapunov function: A function that is non-increasing along all solutions of 
a system of differential equations. Therefore, the corresponding solution 
must move downward through the level sets of the function (i.e., the 
sets where the function is constant). Such a function can be used to 
derive the stability properties at an equilibrium without solving the 
underlying equation. 

separatrix: The kind of solution that begins at a saddle point of a 
planar system of differential equations and tends to another such point as 
time goes on. 


Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chap. 5.3 
Guckenheimer and Holmes, Nonlinear Oscillations , chap. 2.2. 
Hirsch, Smale, and Devaney, Differential Equations , chap. 9.4. 
Roberts, Ordinary Differential Equations, chap. 10.7. 

Strogatz, Nonlinear Dynamics and Chaos , chaps. 4.3 and 6.7. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Duffing, Pendulums. 
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Problems 


1. For a typical solution of the ideal pendulum differential equation, at 
which position is the bob moving fastest? 


2. For a typical solution of the ideal pendulum differential equation, at 
which position is the bob moving slowest? 


3. For the ideal pendulum, use linearization to determine the type of the 
upward equilibrium point. 


4. For the ideal pendulum, use linearization to determine the type of the 
downward equilibrium point. 


5. Use linearization to determine the type of the downward equilibrium 
point for the damped pendulum. 


6. Use linearization to determine the type of the upward equilibrium point 
for the damped pendulum. 


7. Is x' = x 2 + y 2 , y = -2 xy a Hamiltonian system? 
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8. Consider the linear system Y = AY, where 


A = 


v 


c 


b\ 

d J 


Determine conditions on a, b, c, and d that guarantee that this system is 
Hamiltonian. What types of equilibrium points can occur for such a 
linear Hamiltonian system? 


9. Find the Hamiltonian function for the linear system in problem 8. 


Exploration 


Duffing’s equation is a model for a mass-spring system with a cubic 
stiffness term. It is also a model for the motion of a magnetoelastic beam 
mounted between 2 magnets. The second-order equation is given by y" + by 
- y + y 3 . First, assuming the damping constant b = 0, show that this is a 
Hamiltonian system. Second, what happens when b is nonzero? Third, use a 
computer to investigate the behavior of the Duffing equation when a 
periodic forcing term, say Fcos(wf), is introduced. 
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Periodic Forcing and How Chaos Occurs 

Lecture 20 


I n this lecture, we turn to one of the most interesting developments in 
the modern theory of differential equations: the realization that 
solutions of systems of differential equations may behave chaotically. 
Our first example is the periodically forced pendulum. We can force the 
pendulum in one of 2 ways: moving the pendulum up and down 
periodically or forcing the bob alternately in the clockwise and 
counterclockwise directions. 

In the first case, the system of differential equations is given by 
6 ] -v 

v ? = -Z?v - g sin(0) - F cos(tf^) sin(<9). 


Here F measures how far we move the pendulum up and down, and In / co 
is the period of the forcing. 

In the undamped case, we see that the solution curves meander around the 
phase plane in a very complicated fashion. Sometimes the pendulum moves 
around in the clockwise direction; at other times, it moves around in the 
opposite direction. The question is can you predict 5 minutes from now 
which way the pendulum will be swinging? The answer is very definitely 
no; this unpredictability is one of the hallmarks of chaotic behavior. 
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Lecture 20: Periodic Forcing and How Chaos Occurs 


Figure 20.1 



Another ingredient of chaos is what is known as sensitive dependence on 
initial conditions. Basically, this means the following: If we start the 
pendulum out at 2 different but very nearby initial points (0 o ,v o ) and 
(9 X , Vj), then by continuity, the corresponding solutions start out in very 
similar fashion. But before long, the 2 solutions diverge from one another, 
and the corresponding motions of the pendulum are vastly different from 
one another. 

In order to see how we understand this chaos, let’s consider the Lorenz 
system of differential equations. This was the first system of differential 
equations that was shown to behave in a chaotic fashion. Edward Lorenz 
(1917-2008) began his career as a mathematician but turned his attention to 
meteorology while in the army during World War II. In an attempt to 
understand why meteorologists had a hard time predicting the weather, 
Lorenz suggested a simplified model for weather. Basically, Lorenz 
suggested trying to predict the weather on a planet that was surrounded by a 
single fluid particle. This particle is heated from below and cooled from 
above, and like the entire atmosphere on Earth, the particle moves around in 
convection rolls. 
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The Lorenz system of differential equations that describes this motion is 
given by 

x' = a(y - x ) 
y' = Rx-y~xz 
z' = —bz + xy. 


Here a , b , and r are parameters. For simplicity, let us set a = 10 and 
b = 8/3. The remaining parameter R > 0 is called the Rayleigh number. 
We can easily check that there are 3 equilibrium points for this system. One 
is at the origin, and the other 2 are given by (x*,x*,,R-l) and 

(-x*, -x*, R - 1), where x* = (R -1) . When R is relatively small, most 

solutions tend to one of the 2 nonzero equilibrium points that correspond to 
the fluid particle moving periodically in either the clockwise or 
counterclockwise direction. 

But as R increases, suddenly the solutions no longer tend to the equilibria. 
Rather, they tend to circle around these 2 points in a kind of haphazard 
fashion. If we choose 2 nearby initial conditions, just as in the forced 
pendulum case, the corresponding solutions very quickly diverge from one 
another, and we again have sensitive dependence on initial conditions. The 
figures below are 2 views of the same orbit: one projected into the x-z plane 
(left), the other into the y-z plane (right). This type of sensitive dependence 
has been called the butterfly effect. 
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Lecture 20: Periodic Forcing and How Chaos Occurs 


Figure 20.2 



Since our solutions lie in 3-dimensional space, it is difficult to see what is 
actually happening. It appears that the solutions are tending toward a 2- 
dimensional object called the Lorenz template (or the Lorenz mask). Note 
that 2 lobes of the template are joined along a straight line, one lobe in front 
and the other in the rear. 


Figure 20.3 
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Then it appears that solutions wind around this template, looping around the 
2 holes and crossing the central line over and over again. 


Figure 20.4 



Technically, this cannot happen because we would then have different 
solutions that merge into the same solution, contradicting the existence and 
uniqueness theorem (which we discussed earlier). But the fact is that 
solutions do come very close to an object that is similar in shape to this 
template—only there are infinitely many leaves that are bundled closely 
together. This object is known as the Lorenz attractor. 

Most solutions on this template return over and over again to the central line 
as depicted above. So instead of looking for the entire intricate solution of 
the Lorenz system, what we can do is look at how this solution returns over 
and over again to this line denoted by L. So we can think of these solutions 
as being determined by the first return function defined on L. This function 
assigns to each point p on L the next point along the solution through p that 
lies on L. So this gives us the first return function F: L —» L. The 
approximate shape of the graph of F is below. 
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Lecture 20: Periodic Forcing and How Chaos Occurs 


Figure 20.5 



Now the solution starting at p can be tracked by iterating this function. That 
is, given p , we first compute F(p) to find the next return to L. Then we 
compute F(F(p)) to get the second return, and so forth. We denote this 
second iterate of F by F 2 (p). Then F 2 (p) gives the third return, and so forth. 
The collection of points F n (p) is called the orbit of p under F, and 
determining how this orbit behaves gives us a good idea of what is 
happening to the orbits on the Lorenz attractor. 

Amazingly, many simple functions have chaotic behavior when iterated. 
For example, look at what happens to the orbits of 0 and 0.0001 for the 
quadratic function x 2 - 2 or the orbits of .5 and .5001 for another quadratic, 
4x(l -x). 


Important Term 


orbit: In the setting of a difference equation, an orbit is the sequence of 
points x 0 , xi, x 2 , ... that arise by iteration of a function F starting at the seed 
value x 0 . 
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Suggested Reading 


Blanchard, Devaney, and Hall, Differential Equations , chaps. 2.5 and 8.5. 
Guckenheimer and Holmes, Nonlinear Oscillations , chap. 2.3. 

Hirsch, Smale, and Devaney, Differential Equations, chap. 14. 

Strogatz, Nonlinear Dynamics and Chaos , chap. 9. 


Relevant Software 


Blanchard, Devaney, and Hall, DE Tools , Butterfly Effect, Lorenz 
Equations, Pendulum Sensitive Dependence, Pendulums. 


Problems 


1. a. Consider the Duffing oscillator equation with no forcing 
y' = v 
v'=y-y 3 . 

What are the equilibrium points for this system? 

b. Show that this system is Hamiltonian with Hamiltonian function 

H(y, v) = v 2 /2 -y 2 / 2 +//4. 

c. Sketch the level curves of this Hamiltonian function. 

d. Describe the motion of the slender steel beam in the different cases 
that arise. 
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2. Show that for the Lorenz system, there is a straight line solution 
converging to the origin along the z-axis. 


3. Show that the function 

L(x, y, z) = x 2 + 10y 2 + 10z 2 

is a Lyapunov function for the Lorenz system when R < 1. 


4. Given the results of the problem above, what can you say about 
solutions when R < 1 ? 


5. Let V(x, y, z) = Rx 2 + 10y 2 + 10(z - 2 R) 2 . Note that V(x, y,z) = c defines 
an ellipsoid centered at (0, 0, 2 R) in 3-dimensional space. Show that we 
may choose c large enough so that V is decreasing along any solution 
curve that starts outside the ellipsoid given by V(x, y, z) = c. 


6. What does this say about solutions of the Lorenz system that start far 
from the origin? 
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Exploration 


Consider the Rossler system given by 

x ? = —y ~ z 
y' = x+y/4 
z' = 1 + z(x - c ) 

where 0 < c < 7. First find all equilibrium points for this system. Describe 
the bifurcation that occurs at c = 1. Then use a computer to investigate what 
else happens for this system as c increases. Are there any similarities to the 
Lorenz system? 
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Understanding Chaos with Iterated Functions 

Lecture 21 


A s we discussed in the previous lecture, the way mathematicians 
have begun to understand the chaotic behavior that occurs in 
higher-dimensional systems of differential equations is often by 
iterating a first return map. 

Here is another way iteration arises. Recall the limited population growth 
model (a.k.a. the logistic population model) that we saw earlier. This was 
the differential equation;;' = ky(l -y/N) that we solved in multiple different 
ways. Here is a variation on this theme—the discrete logistic population 
model. Suppose x n denotes the population of some species in year n (or, for 
example, generation n). The discrete logistic population model is given by 

x n+ i = kx n {\ - xJN). 

That is, the population given in year n + 1 is just kx n (l - x n /N), where k and 
N are constants that depend on the particular species. Here, N is the 
maximal population rather than the ideal population. For simplicity, let’s 
look at the equation 

x n -\~\ kx n (\ x } ^). 

Here x„ represents the fraction of the maximal population. This is an 
example of a difference equation. To find the population in ensuing years, 
we simply start with an initial population (or seed) x 0 and plug it into the 
function X*) = kx( 1 - x). Out comes xi. Then we do the same with x\ to 
generate x 2 . We do this over and over to produce x 0 , x b x 2 , ..., the so-called 
orbit of the seed x 0 . 

For difference equations, let’s consider more generally x n+ i =f[x n ) for some 
function f We iterate / to produce the orbit of some seed x 0 . For example, 
suppose/^) = x. Then the orbit of 0 is 0, 0, 0,..., so we say that 0 is a fixed 
point for f We see that -1 is an eventually fixed point since f[~l) = 1, 
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which is also fixed. If 0 < |x| < 1, then the orbit of x tends to the fixed point 
at 0. But if |x| > 1, repeated squaring sends the orbit of x off to infinity. So 
we understand the fate of all orbits in this simple case. 

As another example, = x 2 - 1, the points 0 and -1 lie on a cycle of 

period 2 since f[~l) = 0 while f(0) = -1. We also say that 0 and -1 are 
periodic points with prime period 2. More generally, x 0 is a periodic point of 
prime period n if x n = x 0 and n is the smallest positive integer for which this 
is true. 

As with earlier parts of this course, we can visualize orbits in several 
different ways. One way to do this is to plot a time series for the orbit. For 
example, for our cycle of period 2 for j{x) = x 2 - 1, the time series 
representation of the orbit would be the following. 


Figure 21.1 



A second way to visualize orbits is to use graphical analysis. Given the 
function /, we plot its graph and then superimpose the line y = x (the 
diagonal line) on this graph. Given a seed x 0 , we first draw a vertical line 
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from the point (x 0 , x 0 ) on the diagonal to the point (x 0 ,X x o)) = (*o ? * 1 ) on the 
graph. Then we draw a horizontal line back to the diagonal yielding (xi, xi). 
We continue to (x b x 2 ) on the graph and over to (x 2 , x 2 ) on the diagonal. For 
example, below is the graphical analysis representation of the orbit of 0.9 
under f{x) = x 2 (left) and the orbit of-0.7 under f(x) = x 2 - 1 (right). Note 
that the orbit of-0.7 in this second graph tends to the 2-cycle at 0 and -1. 


Figure 21.2 




Unlike the logistic differential equation where all solutions essentially do 
the same thing, for the discrete logistic model, orbits can do many different 
things. We will use the orbit of the critical point for the function^) = kx( 1 
- x) from now on. A critical point is a point where the derivative of the 
function vanishes, so when/(x) = kx{\ - x), we have/'(x) = k — 2Ax, so the 
only critical point is x = 1/2. When 0 < k < 1, the critical point tends to 0, 
which is a fixed point. 
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Figure 21.3 



When 1 < k < 3, graphical analysis shows that the critical orbit now tends to 
a nonzero fixed point. This point is the nonzero root of the equation 
hc( 1 — x)=x, namely (k - l)/k. 


Figure 21.4 
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For k = 3.3, the critical orbit now tends to a 2-cycle, but when k = 3.47 the 
critical orbit tends to a 4-cycle. 

Figure 21.5 




When k = 3.83, the critical orbit tends to a 3-cycle. The time series for this 
critical orbit is below. 


Figure 21.6 
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When k = 4, things go crazy. The orbit of 1/2 is quite simple: It is 1/2, 1, 0, 
0, 0, ..so the orbit is eventually fixed. But any nearby orbit behaves vastly 
differently. The orbit of 0.5001 as a time series is below. 


Figure 21.7 



And here it is using graphical analysis. 


Figure 21.8 
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Lecture 21: Understanding Chaos with Iterated Functions 


We find similar behaviors when we iterate the quadratic function x 2 + c. 
When c = 0, the orbit of 0 is fixed; c = -1, a 2-cycle; c = -1.3, a 4-cycle; 
c = -1.38, an 8-cycle; c = -1.76, a 3-cycle. And again we get chaotic 
behavior when c = - 2. 

Clearly, lots of things are happening when we change the parameter k in the 
logistic function or c in x 2 + c. To see all of this, we plot the orbit diagram 
for this function. In this graphic, we plot the parameter k horizontally with 0 
< k < 4. Along the vertical line above a given &-value, we plot the eventual 
orbit of the critical point 1/2. By “eventual,” we mean that we iterate, say, 
200 times but only display the last 100 points on the orbit. That is, we throw 
away the transient behavior. 


Figure 21.9 
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And here is a magnification of the orbit diagram with k in the interval 
3<k<4. 



difference equation: An equation of the form y n+ 1 = F(y n ). That is, given 
the value y n , we determine the next value y n+ \ by simply plugging y n into the 
function F. Thus, the successive values y n are determined by iterating the 
expression F(y). 

fixed point: A value ofy 0 for which F(y 0 ) =y 0 . Such points are attracting if 
nearby points have orbits that tend to y 0 ; repelling if the orbits tend away 
fromy 0 ; and neutral or indifferent ify 0 is neither attracting or repelling. 

orbit diagram: A picture of the fate of orbits of the critical point for each 
value of a given parameter. 
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Suggested Reading 


Alligood, Sauer, and Yorke, Chaos , chap. 1.5. 

Blanchard, Devaney, and Hall, Differential Equations , chap. 8.1. 
Devaney, A First Course in Chaotic Dynamical Systems , chap. 3. 
Hirsch, Smale, and Devaney, Differential Equations , chaps. 15.1-15.3. 


Relevant Software 


Nonlinear Web, http://math.bu.edu/DYSYS/applets/nonlinear-web.html 

Orbit Diagram for Cx(l-x), http://math.bu.edu/DYSYS/applets/bif- 
dgm/Logistic.html 


Problems 


1. a. Find all fixed points for the function F(x ) = x 2 . 

b. What is the fate of all other orbits for this function? 

2. What is the fate of all orbits of F(x ) = x 2 + 1? 

3. a. What is the fate of all orbits of F(x) = ax where 0 < a < 1? 

b. What is the fate of orbits of this function for other values of al 
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4. Find all fixed points for the logistic family kx( 1 - x) with k > 0. 


5. Find the cycle of period 2 for the function F(x) = -x 3 . Is there a 2-cycle 
for the function F(x) = x 3 ? 


6. Let D(x ) denote the doubling function defined on the interval 0 < x < 1 
and given by 

2x if 0 < x < 1/2 and 

2x — 1 if 1/2 <x < 1. 

Show that the seeds 1/3, 1/7, and 1/15 lie on cycles; and compute 
their periods. 


7. Find all points that lie on cycles of prime period 2, 3, and 4 for the 
doubling function. 


8. Draw the graph of D(x), D 2 (x), and D 3 (x). How many times does the 
graph of D n {x) cross the diagonal? Find all points that lie on cycles of 
period n for the doubling function (not necessarily prime period n). 
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Exploration 


Using the applets below, investigate the orbit diagrams for the functions 
x 2 + c, cx(l - x 2 /3), and csin(x). Do you observe any similarities in these 
diagrams as you zoom in? 


Orbit Diagram for x 2 + c: 

http://math.bu.edu/DYSYS/applets/bif-dgm/Quadratic.html 
Bifurcations of the cubic map: 

http://math.bu.edu/DY S Y S/applets/bif-dgm/Cubic.html 
Bifurcations of the sine map: 

http://math.bu.edu/DYSYS/applets/bif-dgm/Sine.html 
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Periods and Ordering of Iterated Functions 

Lecture 22 


A s we saw when we studied first-order differential equations, there 
are 3 different types of equilibrium points—sinks, sources, and 
nodes. The same is true for iterated functions: There are 3 different 
types of fixed points—attracting, repelling, and neutral (or indifferent). 

A fixed point is attracting if nearby seeds to the fixed point have orbits that 
tend to the fixed point. It is repelling if nearby orbits tend away from the 
fixed point. And it is neutral if neither of the above cases occurs. For 
example, X*) = has an attracting fixed point at x 0 = 0 and a repelling fixed 
point at xq = 1. This is easily seen with graphical analysis. 


Figure 22.1 



Also,X*) = x 2 + 1/4 has a neutral fixed point at x 0 = 1/2 since nearby orbits 
to the left of the fixed point tends toward it, while orbits to the right of the 
fixed point tend away from it. 
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Figure 22.2 



We can use calculus to determine the types of fixed points. Suppose we 
have a fixed point x 0 , soX*o) = *o- If |/'(*o)l < 1, then graphical analysis 
shows that this fixed point is attracting. If [f'(x 0 )| > 1, then x 0 is repelling. If 
f’(x 0 ) = ±1, then we get no information; x 0 could be attracting, repelling, 
or neutral. 

If we have a periodic point x 0 of period n, then we know that f n (x 0 ) = x 0 . 
So this periodic point is attracting or repelling depending on whether 
|(f)'(x 0 )| is less than or greater than 1. For example, for the periodic points 0 
and -1 of period 2 for fix) = x 2 - 1, we have / 2 (x) = x 4 - 2x 2 , so (f 2 )'( 0) = 
(f 2 )'(- 1) = 0. Therefore these periodic points are attracting. Graphical 
analysis of both j{x) and/ 2 (x) shows this as well. 


Figure 22.3 
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Bifurcations also arise when we iterate functions. For example, consider 
f{x) = x 2 + c. There is a bifurcation when c = 1/4 as is seen in the graphs 
when we pass from c > 1/4 to c < 1/4. There are no fixed points when 
c > 1/4; a single (neutral) fixed point when c = 1/4; and 2 fixed points 
(1 attracting, 1 repelling) when c is slightly below 1/4. This type of 
bifurcation is called a saddle-node bifurcation, just as in the case of 
first-order differential equations. 


Figure 22.4 




There is a very different type of bifurcation that we saw many times in the 
orbit diagram. This is the period doubling bifurcation. What happens at this 
bifurcation is an attracting periodic cycle of period n suddenly becomes 
repelling, and at the same time, an attracting cycle of period In branches 
away. For example, consider again f(x) = x 2 + c. When c < 1/4, we always 
have 2 fixed points at 

x + = — ±-. 

2 2 

The fixed point at x+ is always repelling. At x~ we have 
f\x_) = 1 — -s/l — 4c . 


213 










Leture 22: Periods and Ordering of Iterated Functions 


So this fixed point is attracting if-3/4 < c < 1/4 and repelling if c < -3/4. 
We have/'(x_) = -1 if c = -3/4. More importantly, graphical analysis shows 
that a new periodic point of period 2 is bom when c goes below -3/4. 
Equivalently, one can check that the equation/ 2 (x) = x has only 2 real roots 
when c > -3/4 (the 2 fixed points) whereas there are 4 roots when c < -3/4. 
This is the period doubling bifurcation. 





One of the most amazing theorems dealing with iteration of functions on the 
real line is a result proved by Alexander Sharkovsky in the 1960s. This 
theorem says the following. Consider the following ordering of the natural 
numbers: First list all the odds (except 1) in increasing order, then list 2 
times the odds (again except 1) in increasing order, then 4 times the odds, 8 
times the odds, and so on. When you have done this, the only missing 
numbers are the powers of 2, which we then list in decreasing order. Below 
is the Sharkovsky ordering. 

3=>5=>7=>9=>... 

2*3 ^ 2*5 ^ 2*7 2*9 ^... 

4*3 =^> 4*5 =^> 4*7 => ... 
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8»3 => 8*5 => 8*7 => ... 


... 2 3 => 2 2 => 2 => 1 


Sharkovsky’s theorem says that suppose you have a function that is 
continuous on the real line. If this function has a periodic point of prime 
period n , then it must also have a periodic point of prime period k for any 
integer that follows n in the Sharkovsky ordering. In particular, if a 
continuous function has a periodic point of period 3, then it must also have 
a periodic point of all other periods! 

Now go back to the orbit diagram for the logistic function. Toward the 
right, we see a little window that had, at least initially, what appears to be a 
cycle of period 3. By Sharkovsky, there must be infinitely many other 
periodic points in this window. Why do we not see them? 


Figure 22.6. The period-3 window in the orbit diagram. 
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Well, the answer is that only the period-3 cycle is attracting; all the other 
cycles are repelling. The reason for this is another amazing fact about 
iteration of polynomials (and, in fact, what are called analytic functions). If 
you have an attracting or neutral periodic cycle, then the orbit of one of the 
critical points must be attracted to it. This means that for the logistic map 
(where there is only one critical point, x = 1/2), we can have at most one 
attracting cycle. Everything else must be repelling. 


Suggested Reading 


Alligood, Sauer, and Yorke, Chaos , chaps. 1.3-1.4. 

Blanchard, Devaney, and Hall, Differential Equations , chaps. 8.2-8.3. 
Devaney, A First Course in Chaotic Dynamical Systems , chap. 5. 
Hirsch, Smale, and Devaney, Differential Equations, chap. 15.2. 


Relevant Software 


Nonlinear Web, http://math.bu.edu/DYSYS/applets/nonlinear-web.html 

Orbit Diagram for Cx(l-x), http://math.bu.edu/DYSYS/applets/bif-dgm/ 
Logistic.html 


Problems 


1. Find the fixed points for F(x) = x 3 and determine their type. 


2. a. Is the cycle of period 2 for F(x ) = x 2 - 1 given by 0 and -1 
attracting or repelling? 

b. Illustrate this using graphical analysis. 
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3. a. For the function F(x ) = ax , what is the type of the fixed point at 0? 
(This type will depend upon a.) 

b. At which a-values does this family of functions undergo 
a bifurcation? 


4. Determine the values of k for which the fixed points in the unit interval 
for the logistic family kx( 1 - x) with k > 0 are attracting, repelling, 
or neutral. 


5. Determine the type of the cycle of period 2 for the function F(x) = -x 3 . 


6. Let D(x) denote the doubling function given by 2x if 0 < x < 1/2 
or 2x - 1 if 1/2 < x < 1. Show that all of the cycles are 
necessarily repelling. 


7. Determine the lvalue for which the logistic function undergoes a 
period doubling bifurcation at the fixed point. 


8. Describe the bifurcations that the function F(x) = ax + x 3 undergoes 
when a = 1 and a = - 1. 
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Exploration 


Using the applets below, investigate the orbit diagram for the functions kx{ 1 
- x) and x 2 + c. In particular, magnify the successive regions that begin 
where period doubling bifurcations occur as displayed below. Do you see 
any similarities? Explain what you observe. 


Orbit Diagram for Cx(l-x): 

http://math.bu.edu/DY S Y S/applets/bif-dgm/Logistic.html 


Orbit Diagram for x 2 + c: 

http://math.bu.edu/DYSYS/applets/bif-dgm/Quadratic.html 
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Chaotic Itineraries in a Space of All Sequences 

Lecture 23 


T hus far we have encountered chaos in a variety of different settings, 
including Euler’s method, the periodically forced pendulum 
equation, the Lorenz system, and iteration of the logistic map. In 
this lecture, we describe how mathematicians begin to understand 
chaotic behavior. 

To keep things simple, let’s concentrate on the logistic map F(x) = 
kx( 1 - x), where k > 4. Here we choose k > 4 so that the maximum value of 
F is greater than 1. When we use the computer, it appears that all orbits go 
to infinity. Of course that is not the case, as there are clearly 2 fixed points. 
Also, the graph of F n indicates that F n has exactly 2 n fixed points. How are 
these cycles arranged, and what are their periods? 

The graph of F shows that all orbits of points x, with x < 0 or x > 1, go to 
-oo. Also, there is a small open interval A\ surrounding 1/2 and containing 
points that are mapped to points to the right of 1. So these points have orbits 
that escape from the unit interval after 1 iteration and then go to oo. Then, by 
graphical analysis, there is a pair of intervals that we call A 2 that contains 
points that map to A u so these points also escape after 2 iterations. Then 
there are 4 intervals A 3 whose points map to A 2 and hence also escape. 
Next, there is a set of T open intervals that contain points that escape after 
n iterations. 

Let Abe the set of points that never escape from I = [0, 1]. So X= I - u A n . 
We want to understand the behavior of the orbits of all points inX We see 
that this set is similar to the famous Cantor middle thirds set. This set is 
obtained as follows. Start with I. Whenever you see a closed interval, 
remove its open middle third. So we first remove from / the open interval 
(1/3, 2/3). That leaves us with 2 closed intervals, so we remove the open 
middle third of each, meaning (1/9, 2/9) and (7/9, 8/9). Then take this 
process to the limit to obtain the Cantor set. 
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The Cantor set is remarkable in many ways. First, it is a fractal set because 
it is self-similar. Second, even though it looks like it contains very few 
points, it actually contains exactly as many as there are in the entire unit 
interval! And the remaining points are not just endpoints of the removed 
intervals. In fact, most points in the Cantor set are not such endpoints. This 
is a strange set indeed. 

To understand the behavior of orbits in X , we move out of the realm of the 
real line and quadratic functions and into a seemingly much more 
complicated environment, the space of all possible sequences whose entries 
are either Os or Is together with the shift map. Toward this end, we break up 
the unit interval / into 2 subintervals: the left half given by / 0 , lying to the 
left of A i, and the right half given by A, on the other side of A x . (This is 
similar to what we observed as part of the Lorenz attractor.) Then, given 
any point x 0 in the unit interval, we associate an infinite sequence S(x) = (s 0 , 
s i, s 2 , ...) of Os and Is to x 0 via the rule s n =j if F n (x 0 ) lies in the interval Ij. 
The sequence S(x 0 ) is called the itinerary of x 0 . 

For example, the itinerary of 0 is S(0) = (0, 0, 0, ...) since 0 is a fixed point 
that lies in 7 0 . Similarly, S(l) = (1, 0, 0, 0, ...) since F( 1) = 0. Note that if 
S(x 0 ) = (so, s i, s 2 , ...), then S(F(x 0 )) = (s b s 2 , s 3 , ...). That is, to obtain the 
itinerary of F(x 0 ), we just drop the first entry in the itinerary of x 0 . This is 
what we call the shift map a. So we have g(s 0 , s u s 2 , ...) = (s i, s 2 , s 2 ,...). 

Now let £ denote the set of all possible sequences of 0s and Is modulo the 
identifications mentioned above; £ is called the sequence space on 2 
symbols. We say that 2 itineraries are close if they agree on the first n 
digits. Two itineraries are even closer together if they agree on the first n + 
d entries. For example, the itineraries (0, 0, 0, 0, 0, ...) and (0, 0, 0, 1, 1, 1, 
...) are fairly close together, but the sequences (0, 0, 0, 0, 0, ...) and (0, 0, 0, 
0, 0, 0, 0, 0, 0, 1,...) are even closer together. 

While the sequence space looks a little crazy, note that we know a lot about 
the fate of orbits of a. For example, we can immediately write down all of 
the cycles of period n for a. Just take any string of digits of length n and 
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repeat this string over and over infinitely often. The resulting sequence in £ 
then lies on a cycle of period n for a. 

The amazing fact is that the dynamics of F on the unit interval are the same 
as that of a on £. That is, each point x in / has a unique associated itinerary 
and vice versa: Given any itinerary in £, there is a unique point in / with 
that given itinerary. The proof of this is not too difficult and is usually given 
in introductory courses on iterated functions and chaos. 

Because of all this, we now have a very good idea about sequences whose 
orbits under the shift map have very different types of behavior. Consider 
the sequence that begins (0, 1, 0, 0, 0, 1, 1,0, 1, 1, ...); in other words, we 
first list all possible strings of 0s and Is of length 1. Then we list all 4 
possible strings of length 2, then we continue with the 8 strings of length 3, 
the 16 strings of length 4, and so forth. This sequence lies in £ and so 
corresponds to a unique point x 0 in I. But look at what happens to the orbit 
of S(xq) under the shift map. If you shift this orbit a large number of times, 
you can arrange that this orbit comes arbitrarily close to any point in £! As 
a consequence, the orbit of x 0 must come arbitrarily close to any other point 
in X. This is what we call a dense orbit. These are exactly the types of orbits 
that seemed to occupy so much of the orbit diagram, and these are the 
typical orbits we see when we choose a random seed in X. 

Do you see sensitive dependence on initial conditions? Given any x in X, x 
has the associated sequence S(x) = (s 0 , s u s 2 , ...)• But any nearby sequence 
has an itinerary that must differ from that of x at some iteration, so the 2 
orbits eventually reach different intervals, / 0 and I x . Indeed, if we just 
change the tail of the itinerary of x at each stage after the n th , we find a 
nearby orbit whose eventual behavior is vastly different. 

This process can be used to analyze chaos in lots of different settings. For 
example, consider the function on the unit circle that simply doubles the 
angle of a given point, that is, 0 -» 20. We can also write this as the 
complex function F(z) = z. Note that 1/3 —» 2/3 -» 1/3, so 1/3 lies on a 2- 
cycle; 1/7 lies on a 3-cycle; and 1/15 lies on a 4-cycle. We can do symbolic 
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dynamics as before by breaking the circle into 2 sets, I 0 = [0, ti) and I\ = [/r, 
0). Similar procedures to the above show that this function is chaotic on the 
unit circle. We could show that 0 —» 30 is chaotic on the unit circle by 
breaking the circle into 3 arcs—[0, 2^/3), [2;r/3, 4^/3), and [4^/3, 2 7r )—and 
then doing symbolic dynamics on the 3 symbols 0, 1, and 2. 


Important Terms 


itinerary: An infinite string consisting of digits 0 or 1 that tells how a 
particular orbit journeys through a pair of intervals / 0 and 7j under iteration 
of a function. 

shift map: The map on a sequence space that just deletes the first digit of a 
given sequence. 


Suggested Reading 


Alligood, Sauer, and Yorke, Chaos , chaps. 1.5-1.8. 

Devaney, A First Course in Chaotic Dynamical Systems , chaps. 9-10. 
Hirsch, Smale, and Devaney, Differential Equations, chap. 15.5. 


Relevant Software 


Iteration Applet, http://math.bu.edu/DYSYS/applets/Iteration.html 
Nonlinear Web, http://math.bu.edu/DYSYS/applets/nonlinear-web.html 
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Problems 


1. a. In the Cantor middle-thirds set, list the 4 intervals of length 1/27 
that are removed at the third stage of the construction. 

b. What is the total length of all of the removed intervals in the first 3 
stages of this construction? 

c. What is the length of all the intervals that are removed at stage 4? 

d. What is the formula for the length of all of the intervals removed at 
stage nl 

e. Using an infinite series, add up the lengths of all removed intervals 
to determine the length of the Cantor set. 


2. Using symbolic dynamics, describe the itineraries of all points that 
eventually land on 0. How many such points are there? 


3. Which points in / have orbits that eventually land on 0? 


4. Give an example of a different point whose orbit fills the setX densely. 


5. Give an example of a sequence whose orbit under the shift map 
comes arbitrarily close to the sequence (0, 0, 0, ...) but never lands on 
this sequence. 


223 




Lecture 23: Chaotic Itineraries in a Space of All Sequences 


6. Consider the shift map on the sequence space of 3 symbols—say, 0, 1, 
and 2. How many periodic points of (not necessarily prime) period n 
are there? Can you write down a sequence that has a dense orbit? 


Exploration 


Consider the doubling function on the interval [0, 1) defined by D(x) = 2x if 
0 < x < 1/2 and D(x) = 2x - 1 if 1/2 < x < 1. (Here we exclude the point 1 
from the domain.) First use a computer, a calculator, or the Iteration Applet 
at http://math.bu.edu/DYSYS/applets/Iteration.html to compute various 
orbits for the doubling function. What do you see? Is this correct behavior? 
Can you explain what the computer is doing? 

Actually, the doubling function behaves just as chaotically as the function 
4x(l - x) does. Use symbolic dynamics to see this by setting I 0 = [0, 1/2) 
and 7i = [1/2, 1). Incidentally, there is another name for the itinerary 
associated to each x. What is this symbolic representation? 


224 





Conquering Chaos—Mandelbrot and Julia Sets 

Lecture 24 


T he goal of this final lecture is to give a brief glimpse at how 
mathematicians are coming to grips with the compete picture of 
systems with chaotic behavior, and in particular, how these systems 
evolve as parameters vary. 

The simplest example of an iterated map that exhibits a wealth of chaotic 
behavior is the logistic family F(x ) = kx{ 1 - x). The big question is whether 
we understand everything that is happening as the parameter k varies. For 
example, how and when do the various periodic windows arise? How do the 
chaotic regimes arise and change? These questions were only finally 
answered in the mid-1990s, and the way we found the answer to this 
question was to pass to the complex plane and iterate complex rather than 
real functions. 

For historical as well as technical reasons, we will consider the quadratic 
function F{z) = z + c rather than the complex logistic maps. Here z and c 
are both complex numbers. The iteration of z 2 + c arose in modern times 
thanks to the foresight of Benoit Mandelbrot. Mandelbrot used the computer 
to plot what is now called the Mandelbrot set (as well as the Julia sets) for 
these maps. The images he produced around 1980 have had a major impact 
on mathematics. First of all, these fractal images are amazingly intricate and 
beautiful, and as a consequence, they led to quite a bit of curiosity in the 
scientific world. More importantly, these images allowed mathematicians to 
use tools from complex analysis to investigate real dynamics. 

The simplest example of a complex iteration occurs when c = 0 (i.e., 
iteration of the complex function F{z) = z 2 ). It is easily seen that if M < 
then the orbit of z tends to 0, which is an attracting fixed point. If, on the 
other hand, |z| > 1, then the orbit of z tends to oo. Finally, if |z| = 1, then the 
orbit of z stays on the unit circle in the plane. Indeed, the behavior of F on 
this circle is quite chaotic. We certainly see sensitive dependence in any 
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Lecture 24: Conquering Chaos—Mandelbrot and Julia Sets 


neighborhood of a point on the circle since, arbitrarily close by, there are 
points whose orbits go far away, either to 0 or to oo. 

For the complex quadratic function F(z) = z + c, the set of seeds z 0 for 
which the orbit does not tend to infinity is called the filled Julia set. This 
set is named for the French mathematician Gaston Julia, who pioneered the 
study of complex iteration back in the 1920s. It turns out that there are only 
2 different types of filled Julia sets for z + c : Either the filled Julia set is a 
connected set (just one piece), or else it is a Cantor set (infinitely many 
point components, sometimes called fractal dust). 

As for the logistic map, the function F{z ) = z 2 + c has a single critical point, 
this time at 0. Amazingly, it is the orbit of 0 that tells us everything about 
the filled Julia set. For if the orbit of 0 goes to oo, the filled Julia set is a 
Cantor set, but if the orbit of 0 does not tend to oo, the filled Julia set is a 
connected set. This is the result that prompted Mandelbrot to plot the set of 
ovalues in the complex plane for which the orbit of 0 does not escape (and 
so the filled Julia set is connected). This intricate and beautiful image is 
what is known as the Mandelbrot set. 


Figure 24.1 
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It was the geometry of the Mandelbrot set that finally allowed 
mathematicians to understand the real dynamical behavior of x 2 + c, and so 
also the logistic family. There are 2 reasons for this. First, the area of 
mathematics known as complex analysis offers many more mathematical 
tools to study iteration of functions like polynomials. And second, looking 
at objects in the plane gives us many more geometric tools to work with. 

Rather than showing how the Mandelbrot set allows us to understand 
the behavior of x 2 + c for real c- values, let’s look instead at another pattern 
that appears. Attached to the main cardioid are infinitely many smaller 
bulbs. Each of these bulbs contains ovalues for which the corresponding 
orbit of the critical point 0 tends to an attracting cycle of some period. This 
period is the same for all ovalues drawn from a given bulb; this is the 
period of the bulb. So how are these periods arranged as we move around 
the main cardioid? 

We first see that the number of spokes in the antenna attached to the bulb 
gives us the period of the bulb. 


Figure 24.2 



period-3 bulb 


period-5 bulb 


227 





Lecture 24: Conquering Chaos—Mandelbrot and Julia Sets 


We next attach a fraction to each bulb. The denominator is the period of the 
bulb. The numerator can be defined in 3 different ways. (1) The location of 
the smallest spoke in the antenna relative to the principal spoke (the spoke 
that connects to the bulb itself) going in the counterclockwise direction 
gives us the numerator. Above are the 1/3 and 2/5 bulbs. (2) If we plot a 
filled Julia set corresponding to a ovalue from the bulb, the smallest ear 
hanging off the central portion of the set, again counted in the 
counterclockwise direction, also determines the numerator. Below are filled 
Julia sets drawn from the 1/3 and 2/5 bulbs. (3) The rotation number of the 
attracting cycle also gives the numerator. 


Figure 24.3 



It turns out that the bulbs are arranged in the exact order of the rational 
numbers. So we see a lot of connections between the geometry of the 
Mandelbrot set and Julia sets and the dynamical behavior. Proving these 
facts is not that easy, but this nonetheless shows how mathematicians use 
tools from diverse areas of mathematics to understand the behavior of 
complicated systems of differential equations. 
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Important Term 


filled Julia set: The set of all possible seeds whose orbits do not go to 
infinity under iteration of a complex function. 


Suggested Reading 


Devaney, A First Course in Chaotic Dynamical Systems , chap. 17. 
———, The Mandelbrot and Julia Sets. 

Mandelbrot, Fractals and Chaos , chap. 1. 

Peitgen, Jurgens, and Saupe, Chaos and Fractals , chaps. 13-14. 


Relevant Software 


The Quadratic Map, http://math.bu.edu/DYSYS/applets/Quadr.html 

The Mandelbrot Set Iterator, http://math.bu.edu/DYSYS/applets/ 
M-setIteration.html 

Orbit Diagram for x 2 + c, http://math.bu.edu/DYSYS/applets/bif-dgm/ 
Quadratic.html 
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Lecture 24: Conquering Chaos—Mandelbrot and Julia Sets 


Problems 


1. Compute the orbit of i under F(z ) = z and describe its fate. 


2. Compute the orbit of 2 i under F(z) = z and describe its fate. 


3. Compute the orbit of z/2 under F{z) = z and describe its fate. 


4. Which orbits of F(z) = z tend to infinity? 


5. Which orbits of F(z) = z tend to the origin? 


6. For the quadratic function x 2 + c, we saw that 
fixed points on the real line when c > 1/4. What 
complex plane? 


7. What is the fate of the orbit of 0 under z + ft Is the Julia set connected? 


8. Use the computer to investigate what happens to the filled Julia sets for 
ovalues along the path c = -3/4 + iA (where A is a parameter). 


there were no 
happens in the 
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9. Use the computer to look at the portion of the Mandelbrot set that lies 
along the real axis. What object in the Mandelbrot set corresponds to 
the windows in the orbit diagram for x 2 + c? 


10. Discuss the bifurcation that occurs at c = 1/4 along the real axis, but 
now from the complex point of view. 


Exploration 


Consider the bulbs hanging off the period-2 and period-3 bulbs attached 
to the main cardioid of the Mandelbrot set. How are the periods of 
these bulbs arranged? You may use The Mandelbrot Set Iterator software 
located at http://math.bu.edu/DYSYS/applets/M-setIteration.html to find 
these periods. 
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Solutions 


Solutions 


Lecture 1 


1. a. y'(t) = 3t 2 + e'. 

b. y"(t) = 6t + e‘. 

c. = 4 + e‘. 

d. The graph is always positive and increasing. Moreover, e‘ tends to 
0 as t tends to -oo, and to oo as t tends to positive infinity. 



e. t = 0. 
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Solutions are of the form y(t) = t+ C, where C is a constant. 
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Solutions are of the form y(t) = t 2 /2 + C. 













5. All solutions tend to 0 since y' < 0 when y > 0. As before, y(t) = e k \ but 
now k < 0. 


6. The general solution here would by y(t) = t 2 H + constant. 


7. Equilibria occur at y = 1 andy = -1. We havey' > 0 ify > 1 andy < -1, 
whereas y' < 0 if—1 <y< 1. So solutions tend to infinity ify > 1 and to 
-1 ify < 1. 


8. A solution to this differential equation for any constant is 
y(t) = constant. So solutions never move; they stay put. 


9. First write y" = g/m. Then we must have y' = ( g/m)t + A for 
some constant A. And theny(f) = ( \l2){glm)t 2 + At + B, where B is any 
other constant. 


Lecture 2 


1. a. t 3 /3 + t 2 / 2. 
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b. y = 0, 2, and -2 are solutions. 

c. y(t) = 0 when t = 0, 1, and -1. y{t) >0if/ L >lor-l<^<0. 



d. The graph of y{t ) = t 2 - 1 is a parabola opening upward and 
crossing the (horizontal) /-axis at t = 1 and / = - 1. 



e. This graph crosses the horizontal /-axis at / = 0 and t = 2. The 
graph is a parabola opening downward since y(t ) tends to -oo as t 
approaches ±oo. 
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Solutions 


2 . 


3. 


0 


▲ 


4. The only equilibrium point is 0. 


5. There are no equilibria for this differential equation. 
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6. For n even, all solutions with y > 0 tend to oo while all solutions with 
y < 0 tend to 0. When n is odd, again all solutions with y > 0 tend to oo 
while all solutions with y < 0 now tend to -oo. So the answer does 
depend on n. 


7. There are equilibria at -1 and 0. \fy > 0, then y' > 0, so solutions go to 
oo. If — 1 < y < 0, then j/ > 0, so solutions increase to 0. If y < -1, then 
y < 0, so solutions tend to -oo. 


8. The function sin (y) has equilibria at nn for each integer n. Between 0 
and n, sin(y) is positive, so solutions increase to y = n. Between n and 
In, sin(y) is negative, so solutions decrease to y = n. Similar behavior 
occurs in other intervals of length 2 n, solutions always tend to the 
equilibria at y = nn when n is odd and away from the equilibria at 
y = nn when n is even. 


9. What about y' = sin 2 (ny) < ? y r is always positive except at the integers 
where y = 0. 


10. One example would be y' = [y( 1 - y)\ since y f > 0 everywhere except 
y = 0 andy = 1, but of course there are many others, likey' =y 2 ( 1 —y ) 2 . 
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Solutions 


Lecture 3 


1. a. y'(t) = 2 1, so y(t) increases when t > 0. 

b. y(t) decreases when t < 0. 

c. y(t) = (t+ l) 2 so the only root is t = ~\ . 

d. y' = 2t + 2, so this function increases when t>— 1. 

e. 




g. -3sin(3^ + 4). 
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2. The only equilibrium point isy « -1, which is a source. 


3. The only equilibrium point isy = 0, which is a node. 


4. 

v 


o 


v 


5. Equilibria aty = 1 (source) andy = -1 (sink). 


6. The only equilibrium point is at y = 1, and there we have that the 
derivative ofy 3 - 1 is 3y 2 . At 1 we gety' = 3, so 1 is a source. 


7. Technically, the existence and uniqueness theorem does not apply 
when y = 0 since the function [y| is not differentiable there. However, 
we have an equilibrium solution y = 0 there, and all other solutions are 
given by y(t) = Ce f when C > 0 or Ce~ f when C < 0, so we do have 
existence and uniqueness aty = 0. 
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Solutions 


8 . 


The equilibrium points are given by ±A m when A > 0 and 0 when 
A = 0. There are no equilibrium points when A < 0. Since the derivative 
of y - A is 2 y, we have that +A is a source while —A is a sink. 
When A = 0, the equilibrium point at 0 is a node. 


9. The equilibria are at 0 and 1 and the derivative is ^4(1 - 2 y). So for 
y = 0, this point is a source when A > 0 and a sink when A < 0. For 
y = 1, we have a sink when A > 0 and a source when A < 0. When 
A = 0, all points are equilibrium points, so they are all nodes. 


10. This equation has an equilibrium point aty = -1. For y < -1, y' < 0, so 
these solutions tend to -oo. For -1 <y< 1, solutions now increase until 
they hit y = 1, where the slope becomes infinite, so the solutions stop 
there. If y > 1, then y f < 0, so solutions decrease until they hit y = 1, 
when again the slope becomes infinite and solutions stop. 


Lecture 4 


1. The only equilibrium point is aty = A , and this point is a source. 


2. No. 


3. 
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There are no equilibrium points when A is nonzero and infinitely many 
when A = 0; every y- value is an equilibrium point when A = 0. 



4. A bifurcation occurs when A = 0. 


5. The equilibria occur aty = 0 and y = -l/A (as long as A is nonzero). 


6. When A > 0, there is an equilibrium point at y = 0 that is a source. 
When v4 < 0, the equilibrium point aty = 0 is a sink. But when A — 0, all 
points are equilibrium points, so we have a bifurcation at ^4 = 0. 


7. We have equilibria aty = 0 andy = A for each A. When A = 0, there is a 
single equilibrium point at y = 0, a node, with all other solutions 
decreasing. When A > 0, the equilibrium point aty = A is a sink while 0 
is a source. The opposite occurs when A < 0. 


8. Only at those that are nodes. 


9. When B = 0, we have a similar situation to that in problem 7. When B > 
0, we always have a pair of equilibria for each A, given by 

a±ja 2 +4b 

2 

When B < 0, these 2 equilibria only exist if A 2 + 4B > 0 (i.e., for 
B > -^4 2 /4). When B = ~A 2 /4, there is a single equilibrium point at 
A/2. When there are 2 equilibria, the larger one is a sink and the other 
is a source. 
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Solutions 


10. When B = 0, we have a single equilibrium point at y = 0 when A < 0. 
This equilibrium point is a sink. There are 3 equilibria if ^4 > 0: one at y 
= 0, a source; and 2 sinks at y = ±4 1/2 . When B ^ 0, the graph of 
F(y) = B + Ay — y 3 has derivative equal to 0 when 3 y 2 = A. Therefore if 
A < 0, the graph of F(y) is strictly decreasing, so there is always a 
single equilibrium point that is a sink. If A > 0, there are now 2 places 
where F'(y) = 0, at y = ±(4/3) : 1/2 . F(y) has a local minimum at the point 
y = -(4/3) 1/2 and a local maximum at y = +(4/3) 1/2 . For B- values for 
which the value of F at the local minimum is greater than 0, the graph 
of F shows that there is only one equilibrium point, a sink. Similarly, 
there is only one equilibrium point if the value of F at the local 
maximum is less than 0, again a sink. But if the local minimum is less 
than 0 and the local maximum is greater than 0, then there are 3 
equilibrium points. The largest and smallest equilibria are sinks, and 
the other is a source. When either the local maximum or local minimum 
value is 0, then a sink and a source merge to become a node. 


Lecture 5 


1. 

a. 

f/2 + 5 t+C. 


b. 

/ 3 /3 +t 2 + t+C. 


c. 

-r 1 + c. 


d. 

~e * + C. 


e. 

C. 

2. 

t = 

e . 


242 



3. t = ln(4). 

4. y(t) = Ce‘. 



5. y(t) = Ce - 1. 
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Solutions 


6 . 


We have the supposed solution 


y(t) = 


De l 

1 + De* ’ 


First use the quotient rule to compute that 


y\t) = 


De 

(l + De‘f 


Next compute 


yQ-y) 


De' ( De' ' 

\ + De\ l + De‘ y 


De 

(l + De l f 


so we see that this is indeed a solution. 


7. Separating and integrating, we find 
— = f ~^r = f dt-t + c 

y J y J 

so that y(t) = —l/(t + c ). This is not the general solution, since when 
j/(0) = 0, we have 0 = y(0) = ~l/c, or multiplying through by c, we find 
that 0 = -1. Is that true? Obviously not. What’s wrong here? The 
answer is that the solution to the initial value problem y(0) = 0 is just 
the constant function y(t) = 0 (i.e., an equilibrium solution). So the 
solutions -1/(7 + c ) is not quite the general solution; we must also add 
in the solution y{t) = 0 to solve all initial value problems. 
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8. The differential equation here is y' = k(y - 80). This equation is 
separable, so we have 

In | y - 80 |= [ ——— = f kdt -kt + c. 

J y-80 J 

Assuming y(t) > 80, we find by exponentiating both sides that 
y(t) = 80 + e kt + c = 80 + De kt . Since y(0) = 200, we also have 
200 = y(0) = 80 + De° so that D = 120. So our solution so far is 
y(t) = 80 + \20e kt . But we also have y( 1) = 180, so that 180 = 
80 + 120^. Therefore 5/6 = e k or k = ln(5/6). Thus the full solution is 
y(t) = 80 + i20e (ln(5/6)) ' = 80 + 120(5/6)'. 


9. First off, the solution satisfying the initial condition y(0) = 1 is clearly 
the constant solution y(t) = 1 (i.e., one of the 2 equilibrium solutions). 
For the other solutions we can again separate and integrate 


f dy 

J 1 


\ dt 


= t + c. 


To integrate 1/(1 ~y 2 ), note that we can break up this fraction into 


1 _ 1/2 1/2 
1 -y 2 \ + y 1-y 


First assume y(0) = 0. Then our solution lies below the equilibrium 
solution;; = 1, so we may integrate to find 


r( 1/2 

1/2) 

- -+ -- 

+y 

1-y) 
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Solutions 


Exponentiating both sides then yields 



= ke 


or 


1-7. 


= De 2> . 


Solving fory then yields 


Z)e 2t - 
1 + De 2 ' 


which satisfies y(0) = 0 when D = 1. Similar calculations show that this 
expression also solves the initial value problem y(0) = 2 when D = - 3. 


10 . Notice that our solution above includes the equilibrium solution 
y(t) = -1 when D = 0. Also, solving the initial value problem y(0) = A 
shows that A = (D - 1)/(1 + D) so that D = -(A + 1 )/(A - 1), which is 
fine as long as A is not equal to 1. But we know the solution that solve 
the initial value problem y(0) = 1; it is our equilibrium solution y(t) = 1, 
so we must add this solution to the other solution to get the full 
general solution. 
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Lecture 6 


1. a. 3/2. 


b. (yi-yoyih-to). 


c. y = —t+l. 

d. y = 0. 

e. t = 1. 


2. / = 2t, so the slope at t = 1 is 2. So our equation so far \sy = 2t + B. To 
determine B , we know that the point (1,1) lies on this straight line. So 
we have 1=2 + 2?, so B = -\. Thus the equation is y = 2t - 1. 


3. ti = 0.1, andj^i = 1.1. 


4 . t 2 = 0.2, j 2 = 1.21, ^3=0.3,^! = 1.331. 
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Solutions 


5 . 



6. a. Clearly, this solution is just e. 

b. Using a spreadsheet, I calculate thaty(l) = 2.593743 when the step 
size is 0.1. 

c. With step size 0.05,1 findy(l) = 2.653298 and with step size .01 I 
find y{ 1) = 2.704814 (your approximations may differ slightly 
depending on the software you use). Clearly, we are getting better 
approximations of the actual solution. 

d. Using the value ofy(l) = 2.718281, the error when the step size is 
0.1 is 0.124538; when the step size is 0.05, the error is 0.064983; 
when the step size is 0.01, the error is 0.013467. So when we cut 
the step size in half by going from step size 0.1 to 0.05, the error 
decreases by approximately one half. And when we decrease the 
step size by 1/5 when we go from step size 0.05 to 0.01, the error 
also decreases by approximately 1/5. 


7. Any solution that starts above the equilibrium point at y = -1 tends 
toward the value y = 1 where the slope field has infinite slope. Here the 
numerical method “goes crazy,” and we again see chaotic behavior. 
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Lecture 7 


1. Only the origin is an equilibrium point. 


2. All points along the line y = —x are equilibrium points. 


3. a. 



b. It appears that all solutions (except the equilibrium point at the 
origin) move away from the origin along a straight line. 

c. Since this system decouples, solutions are of the form 

x(t ) = kie f 
y{t) = k 2 e. 


249 






Solutions 


4. y = v 

V =y 


5. a. The direction field is always tangent to the circles that are centered 
at the origin and point in the clockwise direction. 

b. One solution is x(/) = cos(t) and y(t) = -sin(t). Another is 
x(7) = sin(t) and y(t) = cos(t). Any constant times each of these is 
also a solution. 


6. a. The vector field is horizontal along the x-axis (pointing to the right 
if 0 < x < 1 and to the left ifx<0orx> 1). The vector field is 
vertical along the y-axis and always points toward the origin. At 
any point off the axes, the x and y directions of the vector field are 
the same as the corresponding directions on the axes. 

b. The equilibria are (0, 0) and (1, 0). 

c. On the x-axis, when x > 0, all solutions tend to 1, and when x < 0, 
all solutions tend to -oo. Meanwhile, on the y-axis all solutions 
tend to 0. So if x 0, the solution tends to (1, 0^ and if x ^ 0, the 
solution tends off to oo to the left. 
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Lecture 8 


1. y"+ 3/ + 2y = 0. 

2. y = v 

v' = -2 y - 3v. 

3. The only equilibrium point lies at the origin. 

4 . -4sin(2^) and -4cos(2^). 

5 . 
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Solutions 


6. We need to solve k x + k 2 = A , ~k x - 2k 2 = B for any given values of A 
and B. Adding these equations yields ~k 2 = A + B, so k 2 = -A - B. Then 
the first equation implies that k x = 2A + B. 

7. a. The characteristic equation is (s + 3)(s + 2), and the general 

solution is k x e~ 3t + k 2 e~ 2t . 

b. We must solve k x + k 2 = 0, ~3k x - lk 2 - 1, which yields k x - -1 
and k 2 = 1. 

c. The graph of y(t) increases at first, then reaches a maximum, and 
then slowly decreases to 0. 

d. Initially the mass moves upward, but then it turns around and 
glides directly back to its rest position. 


Lecture 9 


1. 2e 2t cos(3t) - 3^sin(3 1). 


2. e 2 *(cos(3 1 ) + i sin(30). 


-fr + Vfr 2 - 4k 
2 
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4 . In the clockwise direction. 


5. The characteristic equation is s 2 + 5s + 6 = (s + 3)(s + 2) with roots -2 
and -3, so this system is overdamped. 


6. L’Hopital’s rule says that the limit as t —» oo of tie is the same as the 
limit of the quotient of the derivatives (i.e., lie*). This quotient tends to 
0 as t —» oo. 


7. The roots of the characteristic equation are -1 ± /, so the general 
solution is &ie _ *cos(7) + k 2 e~ t sm(i). 


8. The solution that satisfies y(0) = 0 and /(O) = 1 is y(t) = e^sin(f). The 
derivative is y\t) = -^sin (t) + e~ l cos(0- This derivative vanishes when 
sin(7) = cos (t), so the first positive ^-value where y'(t ) = 0 is when 

t = 7r/4. Then we have y(;r/4) = e _;z/4 sin(;r/4) = 


9. The characteristic equation is s 2 + bs + 1, whose roots are 

-ft±Vft 2 - 4 
2 

So the system is undamped when b = 0, underdamped when 0 < b < 2, 
critically damped when b = 2, and overdamped when b > 2. 
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Solutions 


10 . Clearly y{t) = 1 is a constant solution to the nonhomogeneous equation. 
So the general solution, as in the first-order case, is &iCOs(Y) + 
& 2 sin(f) + 1. That is, the mass just oscillates about a point 1 unit 
removed from the natural equilibrium position. 


Lecture 10 


1. a. y(t) = ke~ f + (l/2)sin(7) - (l/2)cos(f). 

b. All solutions tend to the periodic solution (l/2)sin(/) - (l/2)cos(f). 



d. The solutions are now y(t) = ke f - (l/2)sin(f) - (l/2)cos(f), so 
solutions no longer tend to the periodic solution given by 
-(l/2)sin(7) - (l/2)cos(f) as t increases. 

e. As time tends to -oo, solutions now tend to the periodic solution. 
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2. Make the guess of A cos(f) + B sin(f) to find the solution with A = -1/2 
and B= 1/2. 


3. The general solution is k\ cos (t) + k 2 sin(f) + (\/2)e 


4 . When w = V3, the system is in resonance. 


5. In order for cos (wt) + cos (t) to be periodic, we must find a constant A 
for which 

cos(w(t + A)) + cos(t + A) = cos (wt) + cos(t) 

for all ^-values. In particular, this must be true when t = 0. But 
then we have cos(w^4) + cos(v4) = 2. This implies that we must have 
cos(w^4) = 1 = cos(^4). Therefore we need both wA and A to be integer 
multiples of 2k. So we must have wA = 2nn and A = 2m k for some 
integers n and m, so w = n/m. That is, w must be a rational number. 


6. When w is a rational number, solutions are then periodic in t. 
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Solutions 


Lecture 11 


1. 


"8^ 

v 6 y 


2. a. 

7 ? = 7. 


l-i 

b. 

Only the origin. 

3. 

No. 


4. a. Yes. 
b. No. 


a. 


7 ’: 


1 0 
1 2 


7. 


b. The first solution is x(t) = k x e\ Then we must solve / = 2y + k x e f . 
The equation / = 2 y has general solution y(t) = k 2 e 2t . Therefore for 
the nonhomogeneous equation, we guess Ce l so that C = —k\. So 
our solutions are x(t) = k x e and y(i) = k 2 e 2t - k x e\ 


r k x exp(t) ^ 

~ k j 

rn 

h p ^ 

f°l 

K k 2 exp(2^) - k x exp(7 )j 


-b 

T #V2^ 

v 1 . 
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d. We must be able to solve x(0) = A and y(0) = B for any A and B. 
The first equation gives A = k u and the second then says k 2 ~A=B 
or k 2 = A + B, so this is indeed the general solution. 

e. When k\ = 0, we have a straight line solution along the y-axis 
moving away from the origin. When k 2 = 0, we find a straight line 
solution along the line y = —x again moving away from the origin. 
All other solutions also tend away from the origin. 


Lecture 12 


1. The determinant is -5. 


2. Since the determinant is nonzero, the only equilibrium point is at 
the origin. 


3. The trace is 4, and the determinant is 3. 
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Solutions 


4. a. 



b. It appears there are straight line solutions along the x- andy-axes. 


5. 2and-l. 


6. The eigenvalues are just a and b , since the characteristic equation is 
(a - X)(b -X) = 0. 


7. The eigenvalues are both 0, since this is an upper triangular matrix. 
However, every nonzero vector is an eigenvector since multiplying this 
vector by the matrix yields (0, 0): that is, 0 times the given vector. 

8. The characteristic equation here is A, 2 — (1 + 3yf2)X. So the 

eigenvalues are 0 and 1 + 3V2. The eigenvector corresponding 
to the eigenvalue 0 is given by any nonzero solution of the equation 
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x + 3y = 0, so for instance, the vector (1, -1/3) is one such eigenvector. 
The eigenvector corresponding to the eigenvalue 1+3 is given by 
solving the equation V2 x—y = 0. So one eigenvector is (1, V2 ). 


9. The characteristic equation here is k 2 - 6X + 8 = 0, so the eigenvalues 
are 4 and 2. An eigenvector corresponding to the eigenvalue 4 is any 
nonzero vector of the formy = x, so (1, 1) is an eigenvector for this 
eigenvalue. An eigenvector corresponding to the eigenvalue 2 is any 
nonzero vector satisfying y = -x, so (1, -1) is one such eigenvector. 
Thus the general solution is 


k x e 4 


(\\ 


V 1 / 


+ k^G 


1 

v-1 




y 


The solutions with k x or k 2 equal to 0 are straight line solutions moving 
away from the origin and lying along the lines y = x and y = ~x. All 
other solutions also move away from the origin tangentially to the 
straight line solutions along y = ~x. 


10. The solution satisfying this initial condition is found by solving the 
system of equations 

k x + k 2 = 1 
k x -k 2 = 0 

so k x = k 2 = 1/2 yield this solution. 
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Solutions 


Lecture 13 


1. a. The characteristic equation is ^L 2 +1 = 0. 

b. The eigenvalues are ±i. 

c. One eigenvector associated to the eigenvalue i is 

f 5 ^ 

and an eigenvector associated to ~i is 
f 5 ^ 

2. a. The eigenvalues are 0 and -1. 

b. One eigenvector corresponding to 0 is (1, 0) and corresponding to 
-1 is (-1, 1). 

3. Asa system, y” + by' + ky = 0 may be written 

y = v, V = -Ay - bv or 
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The eigenvalues are roots of X 2 + bX + k = 0 (which is the same 
characteristic equation that we saw for the second order equation), and 
so are given by 

-ft±Vfe 2 - 4k 
2 


4. a. The characteristic equation is X 2 - 2aX + a 2 + b 2 = 0, which has 
roots given by a + ib and a - ib. The eigenvector for a + ib is 
found by solving 

—ibx + by = 0 
—bx - iby = 0, 

so y = ix. One complex eigenvector is therefore (1, /). For the 
eigenvalue a - ib , the equations are 

ibx + by = 0 
~bx + iby — 0, 

soy = —ix. One eigenvector in this case is (1, -/). 

b. For the eigenvalues a ± ib , we have a spiral source if a > 0, a spiral 
sink if a < 0, and a center if a = 0. If a = b = 0, we have real and 
repeated 0 eigenvalues. 
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Solutions 


5. 


a. The characteristic equation is 


X 1 -(l + 3>/2)/l = 0 , 


so the eigenvalues are 0 and 1 + 3%/2 . The eigenvector 
corresponding to 0 is given by x + 3 y = 0 or (3, -1). The 
eigenvector for 1 + 3V2 is given by solving V2 x - y = 0 or 

(1, V2 ). In the phase plane, we have a straight line of equilibrium 
points along the line x + 3y = 0. All other solutions lie on straight 
lines with slope V2 , and these solutions tend directly away from 
the single equilibrium point on this line as time goes forward. 


b. Given the general solution 


K 


f 'X \ 


v- 1 J 


+ k 2 e {M ^ t 


f 


\ 


1 ^ 

V2 y 


we must solve 


3 k\ + k 2 = 1 
~k \ + V 2 k 2 = 0 

so we have V2 k 2 = k x . Then equation 1 implies 3 V 2 k 2 + k 2 = 1 so 

that k 2 = 1/ (1 + 3V2 ) and so k\ = V2 / (1 + 3a/ 2 ). So the solution 
is given by 


V2 

f 3 ^ 

g (l+3^)< 

| 

( 1 ^ 

1 + 3V2 

~b 

1 + 3V2 
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Lecture 14 


1. a. 



b. The eigenvalues are 2 repeated. 

c. Every nonzero vector is an eigenvector. 

d. One of many possible forms for the general solution is 


7 2 , 
k t e 


vO , 


+ 


vV 


e. All nonzero solutions move away from the origin along 
straight lines. 


2. We have repeated eigenvalues given by -1. An eigenvector 
corresponding to -1 is given by solving x + Oy = 0, so one eigenvector 
is (0, 1). Next we solve the equations Ox + Oy = 0 and x + Oy = 1 to find 
the special vector (1,0). Then the general solution is 


k x e 




vV 


+ krjte 


f(\\ 


vV 


+ kr^e 




vOy 
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Solutions 


3. We could solve this using the eigenvalue/eigenvector method as in the 
previous question, but it is much easier to proceed as follows. Our 
equations read x f = 0 and / = x. So we have x(7) = k u and then 
integration yields y(t) = kit + k 2 . 


4. a. We have T = a and D = -a , so the path lies along the line D = -I 
in the ZD-plane. When a = - 4 and a = 0, this line meets the 
repeated eigenvalue parabola T 1 - 4D = 0. When a < -4, we are in 
the real sink region. For -4 < a < 0, we have a spiral sink. And 
when a > 0, we have a saddle point at the origin. 

b. So we have bifurcations at a = -4 (change from real to spiral sink) 
and a = 0 (change from spiral sink to saddle). 


5. Here we have T = a and D = 4, so our path is along the horizontal line 
D = 4. We cross the repeated eigenvalue parabola T 1 - 4D = 0 when 
a 2 = 16, so at a = ± 4. When a < -4, we are in the real sink region. Then 
we cross into the spiral sink region. We cross the D-axis when 
a = T = 0 and then enter the spiral source region, and then finally enter 
the real source region when a = 4. So bifurcations occur when a = ± 4 
and a = 0. 
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Lecture 15 


1. a. The x-nullcline is the x-axis, while the y-nullcline is the y- axis. 


b. 



* 

/ 

\ 


c. The only equilibrium point is the origin. 


d. 





e. Since this system is linear and has real eigenvalues , there is 
only one straight line through the origin containing solutions that 
tend to the origin. 


2. The x-nullclines are given by the y-axis and the line y = —x/3 + 50. The 
y-nullclines are given by the x-axis and the line y = —2x + 100. Most 
solutions that start with x and y nonzero will tend to either the 
equilibrium point at (150, 0) or the one at (0, 100). 
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Solutions 


3. The x-nullcline is given by the y-axis and the line y = 50 - x/2. The 
y-nullclines are given by the x-axis and the line y = 25 - x/6. Most 
solutions will tend to the equilibrium point at (75, 12.5). 


4. The x-nullclines are the lines x = 0 and x = 1, and the y-nullcline is the 
parabola x = y 2 opening to the right. So, there are equilibrium points at 
the origin, (1, 1), and (1, -1). We have x' < 0 when x < 0, so all solutions 
in the left half of the plane tend off to infinity. From the directions of the 
vector field, it appears that (1, 1) is a sink and (1, -1) is a saddle. 


5. a. The x-nullclines are the y-axis and the line y = HA- x/A, and the 
y-nullclines are the x-axis and the line y = 1 + x. There are 
equilibrium points (1, 0), (0, 1), and (0, 0). When A < 1, solutions 
that begin with x and y nonzero tend to (1, 0). When A > 1, 
solutions tend to (0, 1). So a bifurcation occurs when A = 1. When 
A = 1, there is a straight line of equilibria along y = 1 -x. 

b. Now the x-nullclines are the y-axis and the line y = 1 - x, and the 
y-nullclines are the x-axis and y = 1 + Bx. Again the equilibrium 
points are at (1, 0), (0, 1), and (0, 0). Now if B > -1, all solutions 
tend to (0, 1) when the initial conditions are not on the x-axis, 
whereas if B < -1, these solutions tend to (1, 0) (unless the solution 
is on the y-axis). 
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Lecture 16 


1. The partial derivative with respect to x is 2 xy + 3x 2 and with respect to 

y is x 2 + 3 y 2 . 


2. The Jacobian matrix is just the coefficient matrix 

fa b\ 

\ c d , 

3. a. The only equilibrium point is at the origin. 

b. The Jacobian matrix is 

^2x P 

,1 0 / 

At the origin, this matrix becomes 

r 

,1 0 / 

c. The eigenvalues of the Jacobian matrix are ±1, so the origin 
is a saddle. 
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Solutions 


4. 


The only equilibrium point is at the origin, and when linearized, the 
system becomes 


The characteristic equation is X 2 - X +1. The eigenvalues are 


1 ±J-3 
2 ’ 

so the origin is a spiral source. Using a computer, you can see that all 
other solutions spiral toward a periodic solution that surrounds the origin. 


5. One of many possibilities is x' = x 2 , / = y 2 , which has a single 
equilibrium point at the origin. The Jacobian matrix is the zero matrix, 
so both eigenvalues are 0. 

6. The equilibrium points are given by y = 1 and x = ±1. The Jacobian 
matrix is 

' 0 f 

v -2x l y 

At (1, 1) the eigenvalues are (1 d= V— 7) / 2 , so we have a spiral source. 
At (-1, 1) the eigenvalues are 2 and -1, so we have a saddle. 


7. The equilibria are given by (nn , n 12 + nur), where n and m are 
integers. Linearization shows that the equilibrium point is a saddle if n 
and m are both even (or both odd), a sink if n is odd and m is even, and 
a source if n is even and m is odd. The lines x = nn are the x- 
nullclines, so the vector field is tangent to these lines, and solutions 
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remain on them. Similarly, the y-nullclines are the lines y= n 12 + m n , 
and the vector field is again tangent to these lines. These vertical and 
horizontal lines therefore bound squares whose comers are a pair of 
saddles, one sink and one source. In each square, all solutions not on 
the bounding lines tend to the equilibrium that is the sink. 

8 . The equilibrium points are at y = 0 and x = ± J~A , so we clearly have a 
bifurcation at A = 0. Linearization yields the eigenvalues 1 and -2 J~A 
at the equilibrium point ( J~A , 0), so this point is a saddle when A > 0. 

The eigenvalues at the other equilibrium point are 1 and 2 y/~A , so this 
point is a source when A > 0. 


Lecture 17 

1. a. The equilibria are (0, 0) and (1, 1). 
b. The general Jacobian matrix is 


r i — y 

v ~y 1 ~ x > 


At (0, 0) and (1,1) this matrix becomes 



and 


( o -i\ 


c. At (0, 0) we have a source, and at (1, 1) we have a saddle. 
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Solutions 


d. The x-nullclines are x = 0 and y = 1. The y-nullclines are y = 0 and 
x = 1. The different regions are: 



v 



2. a. The coexistence equilibrium point when a = b = 1/2 is given by 
x =y = 800/3. The Jacobian matrix at this point is 

'-2/3 -1/3' 
v —1/3 -2/3/ 

so the eigenvalues are the roots of A, 2 + (4/3)A, + 1/3 =0. These 
eigenvalues are then -1 and -1/3, both of which are negative, so 
this equilibrium point is a sink. 











b. The Jacobian matrix is given by 

'\-xl 200 -ay 1 400 -ax/400 " 

v —by / 400 \-y 1200-bxl 400 y 

At (0, 400), this matrix is 

(1 -a 0^ 


so the eigenvalues are 1 - a and -1. Therefore, we have a sink if 
a > 1 and a saddle if a < 1. At the point (400, 0) we find 
eigenvalues -1 and 1 - b, so this point is a sink if b > 1 and a 
saddle if b < 1. 


3. a. The equilibria are (1, 0), andx = y = 1/2. (Technically, the equation 
for y is not defined if x = 0). 

b. The Jacobian matrix is 

^1 -2x-y -x ^ 
v y 1 / x 2 1 —(2 y/x)/ 

At (1, 0) the eigenvalues are -1 and 1, so we have a saddle. And at 
the other equilibrium point, (1/2, 1/2), we have eigenvalues that are 
complex with negative real part, so this equilibrium is a spiral sink. 

c. The x-nullclines are given by x = 0 and 1 - x = y, and the y- 
nullclines are given by y = 0 and y = x. So the nullclines meet at a 
single point that is not on the axes. In the regions between the 
nullclines, we see that the vector field indicates that solutions 
spiral around the equilibrium point. But that does not tell us the 
complete story, as we could have periodic solutions in this region. 
But the computer shows otherwise. 
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Solutions 


Lecture 18 


1. a. There are no limit cycles since for any point (except the origin), 
the vector field points directly northwest. 

b. All solutions lie along straight lines with slope equal to 1. Since 
the origin is the only equilibrium point, all solutions tend to 
infinity except for those on the line y = x, where x and y are less 
than 0. These solutions tend to the equilibrium point at the origin. 


2. On the circle x 2 +y 2 = 1, the system reduces to 

x' - ~y 
y' = x, 

which is a vector field that is everywhere tangent to the unit circle. 


3. On the unit circle, the vector field is now given by 

x' = 0 
y' = x-y, 

so the vector field points vertically on this circle, which means that the 
unit circle is not a periodic solution. 


4. The second equation says that x = xy/( 1 + x 2 ). Substituting this into the 
first equation yields x - 10 = -4x, or x = 2. Then the second equation 
gives 2 = 2y/5 so thaty = 5. 
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5. a. We have r' = 0 when r= 1, so there is a periodic solution along the 
unit circle. (At the origin, we have an equilibrium point.) Since 
0' = 1, this solution moves counterclockwise around the circle. 
When 0 < r < 1, we have r’ > 0, and when r > 1, r' < 0, so all 
nonzero solutions spiral in to this limit cycle at r = 1 , which is 
therefore stable. 

b. We now have periodic solutions when r 3 - 3r 2 + 2r = 0 and when 
r(r - 2)(r - 1) = 0. That is, r = 1 and 2 are periodic solutions as in 
the previous question. Choosing one point in each of the intervals 
between these circles, we see that r f > 0 when 0 < r < 1; r' < 0 
when 1 < r < 2; and r’ > 0 when r > 2. So r = 1 is a stable limit 
cycle while r = 2 is unstable. 

c. In similar fashion, we have a limit cycle at r = n7r, where n is a 
positive integer. Odd integers yield stable limit cycles; even 
integers yield unstable limit cycles. The rotation is now in the 
clockwise direction. 


6. The roots of ar-r 2 + r 3 = r(r 2 - r + a) are 0 and 
l±Vl-4a 

Besides the equilibrium point at the origin, there are two limit cycles 
when a < 1/4, a single limit cycle when a = 1/4, and no limit cycles 
when a > 1/4. The limit cycle r = r_ is stable; r = r+ is unstable. When 
a > 1/4, r' > 0, so all nonzero solutions spiral out to infinity. 
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Solutions 


Lecture 19 


1. When the pendulum is in the downward position (0 = 0). 

2. When the pendulum is in the upward position (0 = tt). 

3. The Jacobian matrix is 

^ 0 O 


-g cos(6 > ) 0 


, so when 0 = n, we have the matrix 




The eigenvalues for this matrix are +-J~g , so this equilibrium point 
is a saddle. 


4. Now the eigenvalues are ±i-J~g , so linearization does not give any 

information. But we know that the system is Hamiltonian, and the level 
curves surrounding this equilibrium point are ellipses, so this 
equilibrium point is a center. 
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5. Since 0 = lux. the Jacobian matrix is 


^ 0 1 ' 

y~S -b, 

with characteristic equation A 2 + bk + g = 0 and roots 

-b + ^b 2 -4 g 
~2 ' 

If b 2 > 4g, the roots are both real and negative, so we have a real sink. 
If 0 <b 2 < 4g, the roots are complex with negative real part, so we have 
a spiral sink. 


6. The Jacobian matrix now has a +g in the lower left entry, 
so the determinant is now -g < 0. This means the equilibrium point 
is a saddle. 


7. We have dF/dx = 2x = -dG/dy, so this system is Hamiltonian. 


8. Let F(x, y) = ax + by and G(x, y) = cx + dy. This linear system is 
Hamiltonian if dF/dx = a = ~d = -dG/dy. The trace of the 
corresponding matrix is then equal to 0, so we can only have a center or 
a saddle equilibrium point (or repeated 0 eigenvalues). 


9. The Hamiltonian function is given by H(x, y) = by 2 /2 + axy - cx 2 /2 
(plus possibly a constant). 
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Solutions 


Lecture 20 


1. a. (0,0), (1,0), and (-1,0). 


b. Computing the respective partial derivatives 


dH 

dv 


dH 

dy 


■-y-y 


shows that the system is Hamiltonian. 



d. If the beam starts to the right or left of the center of the two 
magnets with small velocity, then it just oscillates back and forth in 
either the left or right region, depending on its starting position. 
But if the beam starts with large velocity, then it will move 
periodically back and forth to the left and right regions. 
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2. The function x(t) = 0, y(t) = 0, and z(t) = Ce 8/3t is such a solution. 


3. We first compute that dL/dt = ~20(x 2 + y 2 - (1 + R)xy) - (160/3)z 2 . 
So we need to show that the term x 2 + / - (1 + R )xy is always 
positive away from the origin. This is certainly true when x = 0. 
Along any other straight line y = Mx , this quantity is equal to 
x 2 (M 2 - (1 + R)M + 1), but the quadratic term M 2 - (1 + R)M + 1 is 
always positive if R < 1. 


4. All solutions must tend to (0, 0, 0) when R< 1. 


5. We compute dV/dt = -20(Rx 2 + y 2 + (8/3)(z 2 - 2 Rz)) = -20(Rx 2 + y 2 + 
(8/3)(z - R) 2 ) = (8/3 )R 2 . Note that the equation Rx 2 + y 2 + (8/3)(z - R) 2 = 
K defines an ellipsoid when K > 0. So if K > (8/3 )R 2 , we have 
dV/dt <0. So we may choose a constant C large enough so that the 
ellipsoid V = C strictly contains the ellipsoid Rx 2 + y 2 + (8/3)(z - R) 2 = 
(8/3 )R 2 in its interior. Then we have dV/dt < 0 on the ellipsoid 
V = C + a for any constant a > 0. 


6. Far enough away on the ellipsoid V= C + a, for example, all solutions 
descend toward the ellipsoid V=C. 


ill 


Solutions 


Lecture 21 


1. a. Oandl. 

b. if M < 1, the orbit of x tends to the fixed point at 0. If |x| > 1, the 
orbit tends to infinity. If x = -1, the orbit lands on 1 after 1 
iteration and so is eventually fixed. 


2. All orbits tend to infinity. 


3. a. All orbits tend to the fixed point at 0. 

b. For each a , 0 is always fixed. So we consider orbits of other 
x-values. If a > 1, all these orbits tend to infinity. If a = 1, all these 
orbits are fixed. If a = -1, all these orbits lie on 2-cycles. 
If-1 < a < 0, all these orbits tend to 0. And if a < — 1, all orbits 
tend to ±oo, alternating between the positive and the negative axis. 


4. The fixed points are 0 and (k - l)/k (which only exists in the unit 
interval if k> 1). 


5. The points 1 and -1 lie on a 2-cycle for -x 3 . There is no 2-cycle for x 3 
since the graph of x 3 is always increasing. 


6 . 
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1/3 -> 2/3 1/3, so 1/3 lies on a 2-cycle. 1/7 2/7 4/7 1/7, so 

1/7 lies on a 3-cycle. 1/15 -► 2/15 ^ 4/15 ^ 8/15 ^ 1/15, so 1/15 lies 
on a 4-cycle. 



7. Only 1/3 and 2/3 have prime period 2. The points 1/7 ... 6/7 have prime 
period 3, and 1/15, ... , 14/15 have prime period 4 (with the exception 
of 5/15 and 10/15, which have prime period 2). 


8. The graph of D crosses the diagonal two times, D 2 four times, D 3 eight 
times, and D n 2 n times. Points of the form k/(2 n - 1), where k is an 
integer and 1 < k < T - 1, have (not necessarily prime) period n. 


Lecture 22 


1. The fixed points are 0 (attracting) and ±1 (repelling). 


2. a. This cycle is attracting. 

b. 
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Solutions 


3. a. If -1 < a < 1, 0 is attracting; if a > 1 or a < -1, 0 is repelling. If 
a = 1, all other orbits are fixed, or if a = -1, all other orbits lie on 
two-cycles. So, in both of these cases, 0 is neutral. 

b. Bifurcations occur at a = 1 and a = — 1 . 


4. The derivative of the logistic function is k - 2kx. So the fixed point at 0 
is attracting for 0 < k < 1 and repelling for k > 1. The other fixed point 
is attracting if 1 < k < 3 and repelling if k > 3. 


5. If F(x) = ~x 3 , then F 2 (x ) = x 9 , and we know 1 and -1 lie on a 2-cycle. 
But the derivative of F 2 (x) is 9x 8 , so this cycle is repelling. 


6. We have D' = 2 at all points, so the derivative of D n is 2 n everywhere, 
so all cycles are repelling. 


7. The derivative of the logistic function at the fixed point (k - 1 )/k is 
equal to -1 when k = 3, so this should be the place where the 
bifurcation occurs. The graph of F 2 then shows the emergence of two 
new fixed points as k increases through 3. 


8. This function always has a fixed point at x = 0, and the derivative at 
this point is equal to a. When a > 1, all other points have orbits that 
tend to infinity. When a < 1, there are two more fixed points at 

±Vl -a . When a = -1, the derivative at 0 becomes -1. For a < -1, 
there is a 2-cycle at a pair of symmetrically located points given by 
+^/-(a +1). 
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Lecture 23 


1. a. (1/27, 2/27), (7/27, 8/27), (19/27, 20/27), and (25/27, 26/27). 

b. The removed lengths are 1/3, 2/9, and 4/27, which add up to 19/27. 

c. At stage 4, the removed intervals have total length 8/81 = 2 3 /3 4 . 

d. At stage n, the removed intervals have length 2 n ~ l /3 n . 

e. The infinite series is 

1 2 2 n ~ l l^f2Y“ 1 1 1 n 

3 9 2 n 3^3 J 31 -% 

so the length of the remaining Cantor set is 0. 


2. Any sequence that ends in all zeros eventually lands on (0, 0, 0, ...) 
under the shift map, so any sequence of the form (s 0 , s u s 2 , ... , 
s n 0000 ...) has this property. 


3. The sequences above correspond to the points that lie at the endpoints 
of the intervals we called A n . 


4. If you just interchange all the blocks of length n in the original 
sequence, you get a new point in the sequence space whose orbit is 
dense. Or you could just take the original sequence and throw in a 0 
before each block. There are countless different ways to find such 
interesting orbits. 
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Solutions 


5. 


For example, (01001000100001000001...). 


6. There are 3” points that are fixed under the n th iterate of the shift map in 
this case; they correspond to all possible blocks of length n , which are 
then repeated infinitely often in the sequence space to get the periodic 
point. As in the case of 2 symbols, if we use the sequence that consists 
of all possible blocks of length 1 in order (i.e., 0, 1, and 2), then all 
possible blocks of length 2 (00, 01, 02, 10, 11, etc.) have a dense orbit. 


Lecture 24 


1. / —» -1 —» 1 —» 1 ..., so this orbit is eventually fixed. 


2. 2/ —> —4 —> 16 —> 256 —» ..., so this orbit tends to infinity. 


3. i/2 —» -1/4 —» 1/256 —» ..., so this orbits goes to zero. 


4. Any point z = x + iy where x 2 +y 2 > 1 (i.e., z lies outside the unit circle) 
has an orbit that tends to infinity. 


5. 
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Any point z inside the unit circle has an orbit that tends to the fixed 
point at the origin. 



6. The fixed points always exist in the complex domain and are given by 
the roots of z - z + c = 0, or 

l + Vl-4 c 
2 

Note that these fixed points are complex when c > 1/4. 


7. The orbit of 0 eventually lands on the 2-cycle given by -1 + i and so 
the filled Julia set is connected. 


8. The filled Julia sets are always a scatter of points when A is non-zero. 
When A = 0, this is the only place where the filled Julia set is 
connected. When \A\ ^ 0 is small, it may appear that the filled Julia set 
consists of a single piece, but changing the number of iterations to be a 
much larger number shows that this is not the case. 


9. Each window in the orbit diagram corresponds to a baby 
Mandelbrot set. 


10. As in question 6, the fixed points still exist when c moves above 1/4. 
But now the filled Julia set immediately shatters into a totally 
disconnected set as soon as c increases above 1/4. 


283 




Types of Differential Equations Cited 


Types of Differential Equations Cited 


FIRST ORDER 


Linear: can be homogenous or nonhomogeneous, autonomous 
or nonautonomous. 



Autonomous 

Nonautonomous 

Homogeneous 

y' + ky = 0 

y' + G(t)y = 0 

N onhomogeneous 

y + ky = 2 

y + G(t)y = 21 


Nonlinear: 

Autonomous: y' = y(l -y) 
Nonautonomous: y f = t 2 +y 2 


SECOND ORDER 


Linear: can be homogenous or non-homogeneous, autonomous 
or nonautonomous. 



Autonomous 

Nonautonomous 

Homogeneous 

y" + by'+ ky = 0 

y' + G(t)y = 0 

Nonhomogeneous 

y" + by' + ky = 2 

/ + G(t)y = 2 
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Nonlinear: 


Autonomous: y" =y 2 
Nonautonomous: y" = t 2 


SYSTEMS 


Linear: Y' = AY = -x' = ax + by 
y = cx + dy 


Nonlinear: x' 

= L 


y = sin(x) 



Autonomous 

Nonautonomous 

Linear 

Y , = AY^\yqxqA is 

a constant matrix 

Y' = AY + (sin(r), cos(O) 

Nonlinear 

Y'=(x 2 , sin(v)) 

Y' = F(Y) = (x 2 + t, sin(v) + cos(f)) 


Both linear and nonlinear systems can be autonomous or 
nonautonomous; however, we did not deal with the nonautonomous 
cases in this course. 
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Using a Spreadsheet to Solve Differential Equations 


Using a Spreadsheet to Solve Differential Equations 




H ere are some tips for creating spreadsheets as introduced in 
Lecture 6: “How Computers Solve Differential Equations.” As of 
this writing, most steps were virtually the same for both Mac and 
PC users, and for Microsoft Excel and OpenOffice; however, do consult 
your software documentation if you have difficulties. 


To enter a formula in a cell: 

First type 6 and then enter the corresponding formula, clicking on the cell 
row/column to indicate the variables. Use * for multiplication. 



A 

B 

C 

1 

10 

5 

3 

2 

=A1*C1 



3 





For a constant reference to a cell, say cell Bl, type “$B$1.” 



A 

B 

c 

1 

10 

5 

3 

2 




3 

=$B$1 
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After hitting Enter: 



A 

B 

C 

1 

10 

5 

3 

2 

30 



3 

5 




This works on a PC with Excel 2007 as well as the spreadsheet program 
available in OpenOffice. 


To fill down a formula: 

Whether in Excel or OpenOffice, grab the lower right comer of the cell or 
cells you want to fill down, then drag it down as far as you wish. (Be 
careful not to move the data from the first cell into the next cell.) In this 
process, cell references will automatically increase, e.g., G3 will change to 
G4 in the next row, to G5 in the second row, and so forth. But constant 
references such as “$F$3” will not change; that is, just insert a “$” before 
the letter, and another “$” before the number, of any cell that you want to 
remain unchanged across calculations to determine more than one cell. 


To insert a chart: 

As demonstrated in Lecture 6, first highlight the data you wish to plot. In 
that lecture, I wanted to plot the t and y values generated by Euler’s method, 
so I highlighted all the entries containing these two values. Then click on 
Insert and choose Chart. Many chart options will be shown. I chose the 
scatter plot for the Euler’s method spreadsheet. In that case, we then had a 
choice of how to connect the dots; select the appropriate type of plot. (Also, 
though not demonstrated in the lecture, you can then modify the size and 
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Using a Spreadsheet to Solve Differential Equations 


color of the lines as well as the scale of the axes by clicking on the object in 
the plot that you wish to modify.) 

Another way to do this in Excel 2007 for PC is to first insert a blank chart 
and then select the data you wish to use. To insert a blank chart, go to the 
Insert tab and select the type of Chart you want. Once you’ve inserted the 
chart, you will be given Design options in the Chart Tools toolbar, and you 
can select the option to Select Data. 


To insert a scrollbar: 

This can be a little complicated, but many people who have seen the 
demonstration become curious about how to do this. In Excel, select View 
—> Toolbars —► Forms. For OpenOffice, select View —► Toolbars —> Form 
Controls. Then choose the scrollbar from this menu. To insert the scrollbar 
into your spreadsheet, highlight the area where you wish to place the 
scrollbar. This could be a horizontal or vertical rectangle. 

For Excel 2007 for PC, you must in addition have the Developer tab 
displayed. To display this tab, click on the circular Microsoft Office button 
in the top left comer of the spreadsheet, which opens a list of options and 
your recent documents, and then click on the Excel Options button at the 
bottom of the window. Once you are in the Excel Options window, select 
Popular from the menu on the left, and check the box for “Show Developer 
tab in the Ribbon.” It should be listed among the “Top options for working 
with Excel.” Click OK. 

Once the Developer tab is displayed in Excel, go to it and click on Insert 
and then select the scrollbar option under Form Controls. To find the 
scrollbar, hover your mouse over the various options and read the 
descriptions that pop up for each icon. Click on the scrollbar icon and then 
click in the spreadsheet where you would like the scrollbar to appear. 
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To resize the scrollbar: 


• For Excel, click on the scrollbar to reveal the points around the 
outline. Place your mouse over the point in the bottom right comer 
of the scrollbar so that an arrow appears. Click on that point and 
drag the edge of the scrollbar in or out until it is the appropriate 
size. The default orientation is vertical, but you can change it to a 
horizontal orientation. To change the orientation from vertical to 
horizontal, follow the same steps as for resizing but drag the 
bottom right edge of the scrollbar up and to the right at the same 
time until its orientation shifts. 

• For OpenOffice, first make sure the Design Mode On/Off button in 
the Form Controls toolbar is in the On mode (the button will 
appear highlighted), and then click on the scrollbar button. Then, 
using your cursor, highlight the cells where you wish to place the 
scrollbar, and it should appear. Following the same steps as above 
for Excel, you can drag it by the comers to resize it. 


To arrange for the appropriate output of the scrollbar: 

For example, in our Euler’s method spreadsheet, I wanted the output 
to be the value of delta t, which moved from .01 to 1.01 in steps of size .01. 
That is, I wanted to have 100 different values of delta t as I manipulated 
the scrollbar. To accomplish this, you need to “format” the scrollbar. 
To do this, while the scrollbar you inserted is highlighted, select 
Format —» Control. 

• For Excel, right-click on the scrollbar and select Format Control. 

• For OpenOffice, right click on the scrollbar and select Control, 
which opens the Scrollbar Properties window. Again, the Design 
Mode On/Off button in the Form Controls toolbar must be in the 
On mode (highlighted). 
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One immediate problem is that a scrollbar only puts out nonnegative 
integers. So you will need to use a little algebra to change the types 
of outputs. 

• The first thing to do is to select a cell link where the output will be 
placed. Since the output is an integer and we want something else, 
I would choose an output cell link that is off the given spreadsheet 
page, such as cell Z25. To assign the output cell, enter that cell 
(say it’s Z25) in the Cell Link field of the Format Control window. 

• For OpenOffice, select the Data tab in the Scrollbar Properties 
window and enter Z25 in the Linked Cell field. 


To choose a minimum and maximum output for the spreadsheet: 

• In our case, we wanted the steps to go from .01 to 1.01 in 100 
different steps. So our minimum value would be 0 and the 
maximum would be 100. If this is not the default in the Minimum 
value and Maximum value fields, enter “0” and “100,” 
respectively. Once you have entered these values, click OK, and 
when you move the scrollbar, you will see the entries in cell Z25 
move from 0 to 100. 

• For OpenOffice, enter these values in the Scroll value min and 
Scroll value max fields in the Data tab of the Scrollbar Properties 
window. For OpenOffice, there is no OK button; simply close the 
Scrollbar Properties window. NOTE: To test the scrollbar 
functionality, make sure the Design Mode On/Off button in the 
Form Controls toolbar is in the Off mode. 
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Since we want the entries to change from .01 to 1.01, we then enter in some 
cell that’s visible, say cell F2, the formula =.01*Z25 + .01 and hit Enter. 



A 

B 

C 

D 

E 

F 

G 

1 








2 






=0.01 *Z25 + 0.01 


3 









After hitting Enter: 



A 

B 

c 

D 

E 

F 

G 

1 








2 






0.01 


3 









Then, as you click down the scrollbar (whether in Excel or OpenOffice), 
you see the entries in cell F2 change from .01 to 1.01 as desired. 

After clicking down once on the scrollbar: 



A 

B 

c 

D 

E 

F 

G 

1 








2 






0.02 


3 
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Using a Spreadsheet to Solve Differential Equations 


After dragging the scrollbar all the way to the bottom: 



292 














Timeline 


■ 


1693 .Isaac Newton’s “fluxions” and 

Gottfried Leibniz’s “calculus of 
differences” offer the first published 
solutions to differential equations. 

1700s.Followers of Leibniz—including 

Daniel Bernoulli, Joseph Lagrange, 
Pierre Laplace, and Leonhard Euler— 
continue the development of 
“differential calculus,” while Newton’s 
“calculus of fluxions” remains more 
influential in England. 


1768 .Euler develops a method for 

approximating the solution to a 
differential equation. 

1835 .William Rowan Hamilton defines the 


special system of differential equations 
in which there is a function that is 
constant along all solutions of the given 
differential equation. 


1841.Carl Jacobi advances the study 

of determinants. 

1858 .Arthur Cayley initiates the use 

of matrices to reduce systems 
of differential equations to 
simpler problems. 
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Timeline 


1866.Jacobi refines Hamiltonian 

functions and pioneers the Jacobian 
matrix for evaluating systems of 
partial derivatives. 

1880s.A geometric approach to nonlinear 

differential equations by Henri Poincare 
initiates the qualitative theory of 
differential equations. 

1890s.Aleksandr Lyapunov defines a function 

that is nonincreasing along all solutions 
of a system of differential equations. 

1901.Poincare’s conjecture for showing the 

existence of a limit cycle is proven by 
Ivar Bendixson. 

1895-1905 .Alternative methods for iterative and 

numerical solutions of differential 
equations are offered by Carl Runge 
and Martin Kutta. 

1918.Gaston Julia and Pierre Fatou pioneer 

iteration theory and lay a foundation for 
work on fractals. 

1931.The differential analyzer, an 

electrically-powered mechanical device 
for solving differential equations, is 
invented by Vannevar Bush and others 
at M.I.T. 
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1942.An M.I.T. differential analyzer with 

2000 vacuum tubes begins to make 
important contributions during WWII. 

1950.The first digital difference analyzer is 

created by Northrop Corporation, a 
leading U.S. aircraft manufacturer. 

1950s.Boris Belousov discovers that certain 

chemical reactions could oscillate rather 
than go to equilibrium. 

1963 .While studying an extremely simplified 

problem in meteorology, Edward 
Lorenz offers the first system of 
differential equations shown to possess 
chaotic behavior. 

1960s.Russian mathematician Anatol 

Zhabotinsky revives Belousov’s work 
on oscillating chemical reactions, with 
their findings beginning to diffuse to 
other countries in 1968. 

1960s.Sharkovsky’s Theorem, describing 

the ordering of cycle periods for 
any iterated function, is published 
in Russian. 

1975 .A paper by Tien-Yien Li and James 

Yorke proves that the presence of a 
periodic point of period 3 implies the 
presence of periodic points of all other 
periods, and therefore chaos. 
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Timeline 


1976.An influential paper published 

in Nature by biologist Robert 
May encourages mathematicians 
to turn from an exclusive emphasis 
on higher-dimensional differential 
equations and give more attention 
to the complicated dynamics in 
low-dimensional difference equations 
such as the logistic difference equation. 

1970s.Benoit Mandelbrot’s work on fractals 

and use of computers extends the work 
of Julia and Fatou on iterated functions, 
drawing worldwide attention to visual 
solutions of real-number problems 
using the complex plane. 

1980s.Eminent mathematicians such as Dennis 

Sullivan, John Milnor, and Curtis 
McMullen make major advances in the 
field of complex dynamics, motivated 
by Mandelbrot's introduction of the set 
that bears his name. 

1998.Swedish mathematician Warwick 

Tucker proves the existence of chaotic 
attractor in the Lorenz equations. 
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Glossary 


beating modes: The type of solutions of periodically forced and undamped 
mass-spring systems that periodically have small oscillations followed by 
large oscillations. 

bifurcation: A major change in the behavior of a solution of a differential 
equation caused by a small change in the equation itself or in the parameters 
that control the equation. Just tweaking the system a little bit causes a major 
change in what occurs. See also Hopf bifurcation, pitchfork bifurcation, 
and saddle-node bifurcation. 

bifurcation diagram (bifurcation plane): A picture that contains all 
possible phase lines for a first-order differential equation, one for each 
possible value of the parameter on which the differential equation depends. 
The bifurcation diagram, which plots a changing parameter horizontally and 
the y value vertically, is similar to a parameter plane, except that a 
bifurcation diagram includes the dynamical behavior (the phase lines), 
while a parameter plane does not. 

carrying capacity: In the limited population growth population model, 
this is the population for which any larger population will necessarily 
decrease, while any smaller population will necessarily increase. It is the 
ideal population. 

center: An equilibrium point for a system of differential equations for 
which all nearby solutions are periodic. 

characteristic equation: A polynomial equation (linear, quadratic, cubic, 
etc.) whose roots specify the kinds of exponential solutions (all of the form 
e At ) that arise for a given linear differential equation. For a second-order 
linear differential equation (which can be written as a 2-dimensional linear 
system of differential equations), the characteristic equation is a quadratic 
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Glossary 


equation of the form 2 2 - T 2 + D, where T is the trace of the matrix and D 
is the determinant of that matrix. 

critically damped mass-spring: A mass-spring system for which the 
damping force is just at the point where the mass returns to its rest position 
without oscillation. However, any less of a damping force will allow the 
mass to oscillate. 

damping constant: A parameter that measures the air resistance (or fluid 
resistance) that affects the behavior of the mass-spring system. Contrasts 
with the spring constant. 

determinant: The determinant of a 2-by-2 matrix is given by ad - be , 
where a and d are the terms on the diagonal of the matrix, while b and c are 
the off-diagonal terms. The determinant of a matrix A (det4, for short) tells 
us when the product of matrix A with vector Y to give zero (A Y = 0) has 
only one solution (namely the 0 vector, which occurs when the determinant 
is non-zero) or infinitely many solutions (which occurs when the 
determinant equals 0). 

difference equation: An equation of the form y n+ \ = F(y n ). That is, given 
the value y m we determine the next value y n+ \ by simply plugging y n into the 
function F. Thus, the successive values y n are determined by iterating the 
expression F(y). 

direction field: This is the vector field each of whose vectors is scaled 
down to be a given (small) length. We use the direction field instead of the 
vector field because the vectors in the vector field are often large and 
overlap each other, making the corresponding solutions difficult to 
visualize. Those solutions and the direction field appear within the phase 
plane. See also vector field. 

eigenvalue: A real or complex number usually represented by X (lambda) 
for which a vector Y times matrix A yields non-zero vector XY. In general, 
an n x n matrix will have n eigenvalues. Such values (which have also been 
called proper values and characteristic values) are roots of the 
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corresponding characteristic equation. In the special case of a triangular 
matrix, the eigenvalues can be read directly from the diagonal, but for other 
matrices, eigenvalues are computed by subtracting X from values on the 
diagonal, setting the determinant of that resulting matrix equal to zero, and 
solving that equation. 

eigenvector: Given a matrix A, an eigenvector is a non-zero vector that, 
when multiplied by A, yields a single number X (lambda) times that vector: 
AY= XY. The number X is the corresponding eigenvalue. So, when X is real, 
AY scales the vector Y by a factor of lambda so that AY stretches or 
contracts vector Y (or ~Y if lambda is negative) without departing from the 
line that contains vector Y. But X may also be complex, and in that case, the 
eigenvectors may also be complex vectors. 

equilibrium solution: A constant solution of a differential equation. 

Euler’s formula: This incredible formula provides an interesting 
connection between exponential and trigonometric functions, namely: The 
exponential of the imaginary number (i-t) is just the sum of the 
trigonometric functions cos(f) and i sin(t). So e {lf) = cos (t) + i sin(t). 

Euler’s method: This is a recursive procedure to generate an 
approximation of a solution of a differential equation. In the first-order case, 
basically this method involves stepping along small pieces of the slope field 
to generate a “piecewise linear” approximation to the actual solution. 

existence and uniqueness theorem: This theorem says that, if the right- 
hand side of the differential equation is nice (basically, it is continuously 
differentiable in all variables), then we know we have a solution to that 
equation that passes through any given initial value, and, more importantly, 
that solution is unique. We cannot have two solutions that pass through the 
same point. This theorem also holds for systems of differential equations 
whenever the right side for all of those equations is nice. 

filled Julia set: The set of all possible seeds whose orbits do not go to 
infinity under iteration of a complex function. 
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first derivative test for equilibrium points: This test uses calculus to 
determine whether a given equilibrium point is a sink, source, or node. 
Basically, the sign of the derivative of the right-hand side of the differential 
equation makes this specification; if it is positive, we have a source; 
negative, a sink; and zero, we get no information. 

fixed point: A value of y 0 for which F(y 0 ) =y 0 . Such points are attracting if 
nearby points have orbits that tend to y 0 ; repelling if the orbits tend away 
fromy 0 ; and neutral or indifferent ify 0 is neither attracting or repelling. 

general solution: A collection of solutions of a differential equation from 
which one can then generate a solution to any given initial condition. 

Hopf bifurcation: A kind of bifurcation for which an equilibrium changes 
from a sink to a source (or vice versa) and, meanwhile, a periodic solution 
is bom. 

initial value problem: A differential equation with a collection of special 
values for the missing function such as its initial position or initial velocity. 

itinerary: An infinite string consisting of digits 0 or 1 that tells how a 
particular orbit journeys through a pair of intervals 7 0 and I x under iteration 
of a function. 

Jacobian matrix: A matrix made up of the various partial derivatives of the 
right-hand side of a system evaluated at an equilibrium point. The 
corresponding linear system given by this matrix often has solutions that 
resemble the solutions of the nonlinear system, at least close to the 
equilibrium point. 

limit cycle: A periodic solution of a nonlinear system of differential 
equations for which no nearby solutions are also periodic. Compared with 
linear systems, where a system that has one periodic solution will always 
have infinitely many additional periodic solutions, the limit cycle of a 
nonlinear system is isolated. 
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limited population growth model with harvesting: This is the same as the 
limited population growth model except we now assume that a portion of 
the population is being harvested. This rate of harvesting can be either 
constant or periodic in time. 

Lyapunov function: A function that is non-increasing along all solutions 
of a system of differential equations. Therefore, the corresponding 
solution must move downward through the level sets of the function (i.e., 
the sets where the function is constant). Such a function can be used to 
derive the stability properties at an equilibrium without solving the 
underlying equation. 

mass-spring system: Hang a spring on a nail and attach a mass to it. 
Push or pull the spring and let it go. The mass-spring system is a 
differential equation whose solution specifies the motion of the mass as 
time goes on. The differential equation depends on two parameters, the 
spring constant and the damping constant. A mass-spring system is also 
called a harmonic oscillator. 

Newton’s law of cooling: This is a first-order ODE that specifies how a 
heated object cools down over time in an environment where the ambient 
temperature is constant. 

node: An equilibrium solution of a first-order differential equation that has 
the property that it is neither a sink nor a source. 

orbit: In the setting of a difference equation, an orbit is the sequence 
of points xq, xi, x 2 , ... that arise by iteration of a function F starting at the 
seed value x 0 . 

orbit diagram: A picture of the fate of orbits of the critical point for each 
value of a given parameter. 

ordinary differential equation (ODE): A differential equation that 
depends on the derivatives of the missing functions. If the equation depends 
on the partial derivatives of the missing functions, then that is a partial 
differential equation (PDE). 
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phase line: A pictorial representation of a particle moving along a line that 
represents the motion of a solution of an autonomous first-order differential 
equation as it varies in time. Like the slope field, the phase line shows 
whether equilibrium solutions are sinks or sources, but it does so in a 
simpler way that lacks information about how quickly solutions are 
increasing or decreasing. The phase line is primarily a teaching tool to 
prepare students to make use of phase planes and higher-dimensional 
phase spaces. 

phase plane: A picture in the plane of a collection of solutions of a system 
of two first-order differential equations: x f = F(x, y) and y' = G(x, y). Here 
each solution is a parametrized curve, (x(t), y(t)) or (y(t), v(f)). The term 
“phase plane” is a holdover from earlier times when the state variables were 
referred to as phase variables. 

pitchfork bifurcation: In this bifurcation, varying a parameter causes a 
single equilibrium to give birth to two additional equilibrium points, while 
the equilibrium point itself changes from a source to a sink or from a sink to 
a source. 

predator-prey system: This is a pair of differential equations that models 
the population growth and decline of a pair of species, one of whom is the 
predator, whose population only survives if the population of the other 
species, the prey, is sufficiently large. 

resonance: The kind of solutions of periodically forced and undamped 
mass-spring systems that have larger and larger amplitudes as time goes on. 
This occurs when the natural frequency of the system is the same as the 
forcing frequency. 

saddle-node bifurcation: In an ODE, this is a bifurcation at which a single 
equilibrium point suddenly appears and then immediately breaks into two 
separate equilibria. In a difference equation, a fixed or periodic point 
undergoes the same change. A saddle-node bifurcation is also referred to as 
a tangent bifurcation. 
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saddle point: An equilibrium point that features one curve of solutions 
moving away from it and one other curve of solutions tending toward it. 

separation of variables: This is a method for finding explicit solutions of 
certain first-order differential equations, namely those for which the 
dependent variables (y) and the independent variables (t) may be separated 
from each other on different sides of the equation. 

separatrix: The kind of solution that begins at a saddle point of a 
planar system of differential equations and tends to another such point as 
time goes on. 

shift map: The map on a sequence space that just deletes the first digit of a 
given sequence. 

sink: An equilibrium solution of a differential equation that has the property 
that all nearby solutions tend toward this solution. 

solution curve (or graph): A graphical representation of a solution to the 
differential equation. This could be a graph of a function y(t) or a 
parametrized curve in the plane of the form (x(£), 

source: An equilibrium solution of a differential equation that has the 
property that all nearby solutions tend away from this solution. 

spiral sink: An equilibrium solution of a system of differential equations 
for which all nearby solutions spiral in toward it. 

spiral source: An equilibrium solution of a system of differential equations 
for which all nearby solutions spiral away from it. 

spring constant: A parameter in the mass-spring system that measures how 
strongly the spring pulls the mass. Contrasts with the damping constant. 

steady-state solution: A periodic solution to which all solutions of a 
periodically forced and damped mass-spring system tend. 
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Taylor series: Method of expanding a function into an infinite series of 
increasingly higher-order derivatives of that function. This is used, for 
example, when we approximate a differential equation via the technique of 
linearization. 

trace: The sum of the diagonal terms of a matrix from upper left to 
lower right. 

vector field: A collection of vectors in the plane (or higher dimensions) 
given by the right-hand side of the system of differential equations. Any 
solution curve for the system has tangent vectors that are given by the 
vector field. These tangent vectors (and, even more so, the scaled-down 
vectors of a corresponding direction field) are the higher-dimensional 
analogue of slope lines in a slope field. See also direction field. 
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